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ONS-DIMENSICNAI  HYDRODYNAMIC  CODE  FOR  NDCLEAR-fiXPLOSICK  CALCULATIONS 


A  number  of  requests  ha  vs  been  received  by  this  Laboratory  for 
complete  details  of  the  one -dimensional  hydrodynamic  code  used  ior 
the  nuclear-explosion  calculations  of  references  3  and  4.  This 
report  has  been  prepared  to  satisfy  these  requests.  The  FORTRAN  IV 
listings  given  in  the  appendices  should  be  sufficient  to  enable  the 
reader  to  duplicate  this  code,  if  he  so  desires. 

This  version  of  the  WUNDT  code  has  been  written  explicitly  for  publica^ 
tion;  it  is  not  the  version  actually  used  for  references  3  and  4.  A 
number  of  errors  may  have  escaped  detection  of  the  authors  in  the 
print-outs  of  this  report.  The  authors  welcome  notification  of  any 
errors  found. 

Support  for  development  of  the  WUNDY  (ref.  1)  code  for  nuclear  explo¬ 
sions  was  provided  by  the  Defense  Atomic  Support  Agency  through 
Nuclear  Weapons  Research  Subtask  01.003  (Task  N0L-181) . 


R.  E.  ODENING 
Captain,  USN 
Commander 


C.J.  ARONSON 
By  direction 


ii 


HOLTS  62-168 


comers 

Pag© 

1.  IFBRODOCTIQN. . . . . . . . .  1 

2o  DESCRIPTION  OF  TBS  WUNDT  CODE. . . . . . 

2.1  General  Properties . . . . . 

2.2  The  Finite  Difference  Equations . . . 

2.3  Initialisation......... . . 

2.4  Reaoniag . . . . . . 

2.5  Suss&ry  Routines . . . . . . . 

2.6  Equations  of  State . . . 

3.  USING  THE  CODE. . . . . . . 

3.1  Initial  Conditions  for  Nuclear  Explosion  Problems 

3.2  Other  Problems . . . . . . 

4.  REFERENCES . . . . . . .  7 

APPENDIX  A  -  Examples  A-l 

APPENDIX  B  -  Listings  of  the  WJHDY  Hydrocode  B-l 


iii 


O  £*-  -S'-  Vrf  W  M  »-*  M 


KOLTR  62-168 


1.  INTRODUCTION 

This  s«port  describes  the  one-disasneioaal  hydrodynamic  cods,  WUHDY, 
that  is  being  used  at  the  Haval  Ordnance  Laboratory  for  nuclear-explosion 
calculations. 

Tbs  first  FSSitTRAN  version  of  WUNDY,  written  by  W.  Walker  (rsf.  l), 
was  based  on  the  SO-CODE  of  tbs  Ur1 .  rsity  of  California  Radiation 
Laboratory  (ref.  2).  Several  vex  *-;-aa  of  WUNDY  now  exist  at  this 
Laboratory.  The  one  described  in  this  report  has  been  used  in  the 
calculation  of  some  of  th®  hydrodynamic  aspects  of  nuclear  explosions  in 
air  (ref.  3,  4). 

2.  DESCRIPTION  OF  THE  WUNDT  CODE 

2.1  General  Properties;  A  ona-disensiomi  hydrodynamic  code  is 
basically  quite  aia^jle.  It  is  the  refinements  *hat  sake  it  esea  lengthy. 
Many  options  are  usually  included;  for  exax^le,  pic  no,  cylindrical  and 
spherical  symmetries  are  usually  included  in  the  sea®  code.  Features  life® 
a  variable,  internally  calculated  tisae  step  are  desirable.  Special 
computations  are  often  included,  such  as  the  x-ray  deposition  of  FUFF 
(refs.  5,  6,  7,  8).  Equations  of  state  can  become  quite  long.  The  input 
and  output  routines  usually  occupy  a  good  deal  of  apace. 

WUNDY  as  given  here  is  in  cgs  units.  It  handles  only  pure  hydro¬ 
dynamics  (no  radiative  transfer,  conduction,  etc.)  in  plan©,  cylindrical 
or  spherical  geometry.  The  geometry  of  the  code  is  collected  into  a 
single  accessible  place  where  it  can  be  modified  for  calculations 
requiring  area  to  be  some  function  of  radius  other  than  A  -  cr0  with 
n  =  0,  1,  or  2.  Shocks  are  handled  by  the  artificial  viscosity  method. 
WUHDT  has  300  zones,  divided  among  as  ma ay  as  30  different  regions  or 
mat© rials.  Tbs  version  described  in  this  report  is  written  for  clarity 
rather  than  for  maximum  computing  efficiency.  For  example,  constants  of 
proportionality  (e.g.,  pi)  are  retained  even  though  they  would  eventually 
cancel  out.  Thus,  in  spherical  geometry,  the  "Bass'8  of  a  son®  is  always 
the  actual  number  of  grams  of  material  in  that  sons. 

A  model  of  the  earth's  atmosphere  (ref.  9)  is  included  as  a  subroutine 
and  can  be  used  to  furnish  initial  ambient  atmospheric  conditions. 
Calculations  can  be  performed  for  shockwave  propagation  in  an  arbitrary 
direction  from  an  explosion  at  an  arbitrary  altitude  of  burst  (e.g.,  the 
shockwave  propagation  down  to  the  ground  from  a  high  altitude  explosion 
can  be  calculated).  The  component  of  the  earth's  gravity  along  the 
calculated  ray  is  included  in  the  equations  of  motion. 

2.2  The  Finite  Difference,  Eamtlpjass  The  problem  is  advanced  from 
tics  tn  to  time  tn*l  by  the  usual  finite  difference  equations  (refs.  10, 
11,  12,  13): 
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where  x  *  distance  from  origin 
u  »  velocity  of  interface 
v  *  specific  volume  of  zone 
p  *  pressure  in  zone 
q  "  artificial  viscosity  in  zone 
P  -  p  +  q 

e  “  internal  energy  in  zone 

Y  *  equation  of  state  constant  *  1.4  for  ideal  air 
g  *  acceleration  of  gravity  component 

*  area  of  interface  j  at  time  n,®  4n(x*j)3  for  sphere 

^yi/2  *  VOiume  0f  2°ne  J-l/2  *  |fT^(Xj)5*’(XJ_x)a j  f>°r  8phere 
Itj  ®  timestep  *  l/2  (at + 


m  =  mass  o  1'  zone 


These  equations  are  discussed  in  the  references;  this  material  will 
not  be  repeated  here.  The  equations  for  p  and  e  are  for  a  gansna- 
law  gas,  e  =  pv/(Y-l). 
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2.3  Initialization:  A  problem  is  initialized  when  each  interface 
has  been  assigned  a  position  and  Telocity,  each  zone  between  interfaces 
has  been  assigned  a  mass  and  internal  energy,  and  all  of  the  control 
variables  have  been  specified.  Most  problems  are  completely  initialized 
by  a  deck  of  about  18  cards.  The  first  two  cards  contain  lettering  to  be 
printed  on  the  output.  The  next  four  cards  contain  the  control  variables 
(printing  frequency,  type  of  geometry,  etc.).  These  are  followed  by  a 
group  of  cards  (one  for  eadh  region)  describing  each  region  of  the  problem 
(a  region  is  usually  a  group  of  zones  of  the  same  material) .  The 
dimensions,  velocities,  densities,  number  of  zones  in  each  region,  stc. 
are  given  on  these  cards.  The  normal  input  deck  ends  with  ten  cards 
giving  some  additional  data  (direction  of  calculated  ray,  altitude  of 
burst ,  etc . ) . 

In  some  cases  one  wishes  to  give  initial  data  for  interfaces  or 
zones  individually  rather  by  regions.  For  example,  in  a  moving-piston 
problem,  the  piston  would  be  introduced  by  giving  the  first  interface  a 
velocity.  This  is  done  by  additional  cards  following  the  normal  input 
deck. 

Additional  input  can  also  be  read  from  tape.  Special  input 
calculations  (e.g.,  in  high  explosive  problems,  the  initial  conditions 
in  the  detonation  wave)  can  be  made  by  other  programs  and  written  on 
tape  that  V/UNDI  can  read.  WUKBY  can  also  write  its  own  tape  from  any 
record  of  which  the  problem  can  be  resumed  later. 

2.4  Re  zoning:  The  computer  time  required  for  a  problem  is  about 
proportional  to  the  number  of  zones  that  are  being  calculated  and  about 
inversely  proportional  to  the  width  of  the  smallest  zone  in  the  problem. 
Four  types  of  rezoning  have  been  used  in  reducing  computing  times.  A 
general  rezoning  routine  that  contains  all  four  of  these  methods  as 
special  cases  could  probably  be  devised. 

The  first  type  of  re  zoning  becomes  alerted  when  the  timestep 
begins  to  decrease  appreciably.  In  some  problems  this  decrease  is  due 
to  excessive  compression  of  a  zone  that  does  not  lie  in  the  region  of 
interest.  This  rezoning  method  (REZ1)  seeks  out  the  zone  causing  the 
small  timestep  and  removes  it  by  combining  it  with  an  adjacent  zone. 

A  second  type  of  re zoning  (R2Z2)  waits  until  a  disturbance 
(usually  the  main  shockwave)  reaches  the  outer  boundary  of  the  problem 
and  then  moves  this  boundary  out  by  a  factor  of  2  in  distance.  It  then 
combines  pairs  of  adjacent  zones  to  release  enough  zones  to  fill  the  new 
extra  region  of  space  out  to  the  new  boundary.  This  allows  the  same 
number  of  zones  to  cover  an  indefinitely  expanding  space,  at  the  cost 
of  resolution  lost  by  making  the  zones  larger.  This  is  the  type  of 
rezoning  used  in  the  calculations  of  references  3  and  4. 
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A  third  type  of  re  zoning  (REZ3)  keeps  down  the  total  number  of 
zones  required  for  a  problem  by  having  fairly  large  zones  everywhere 
except  in  the  region  of  interest  (usually  the  main  shockwave) ,  The 
shockwave  is  kept  supplied  with  fine  zones  by  subdividing  large  ambient 
zones  as  the  shock  approaches  and  then  recombining  them  into  large  zones 
when  the  shock  has  passed.  This  scheme  is  now  in  use  on  nearly  every 
problem. 


A  fourth  type  of  rezcning  (KEZ4)  replaces  the  zones  of  the 
first  (fireball)  region  by  a  smaller  number  of  zones.  In  some  problems 
the  first  region  holds  down  the  timestep  needlassly  long  after  the 
shockwave  has  moved  away.  In  these  cases  the  problem  can  be  made  to 
run  faster  by  a  factor  of  two  or  three  by  replacing  the  spent  fireball 
by  just  a  few  large  zones. 

2.5  Summary  Routines:  Two  subroutines  have  been  included  to 
provide  data  in  addition  to  the  normal  output.  Subroutine  0UT2  prints 
tables  of  the  shock  location,  pressure,  etc.  as  a  function  of  time. 

Thi3  output  is  useful  in  preparing  shock  front  pressure  versus  distance 
curves  and  shock  front  radius  versus  time  curves.  Subroutine  0UT3 
prints  tables  of  the  pressure,  density,  particle  velocity,  dynamic 
pressure,  etc.  versus  time  at  fixed  distances  (compared  to  the  normal 
output,  which  gives  quantities  versus  distance  at  fixed  times) .  This 
output  is  useful  for  predicting  the  pressure-time  records,  etc.  at 
fixed  gage  locations.  An  option  is  included  for  collecting  data  at 
points  moving  in  straight-line  trajectories  near  the  explosion. 

2.6  Equations  of  State :  Two  equations  of  state  are  included  in 
the  listings  of  Appendix  B.  One  (EQSl)  is  for  an  ideal  gas  of  constant 
gamma.  The  other  (EQS2)  is  an  approximate  equation  of  state  for  real 
air,  based  on  reference  1 4.  These  can  be  used  as  patterns  for  preparing 
equations  of  state  for  other  materials.  The  essential  things  to  be 
calculated  are  the  pressure,  internal  energy,  and  an  approximate  sound 
speed  (needed  in  calculating  the  timestep) .  The  temperature  and  gamra 
may  be  calculated  if  desired  but  are  not  needed  elsewhere  in  the 
program. 

3.  USING  THE  CQIE 

3.1  Initial  Conditions  for  Nuclear  Explosion  Problems;  One 
approximation  to  a  nuclear  explosion  that  has  been  found  useful  (refs. 

3  and  4)  is  a  uniform  sphere  of  highly  heated  air  that  is  allowed  to 
expand  hydrodynamically  into  the  surrounding  ambient  air. 
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The  energy  of  the  explosion  (adjusted,  if  desired,  for  early 
losses  by  thermal  radiation)  is  assumed  to  be  uniformly  distributed 
within  a  sphere  whose  radius  in  meters  can  be  taken  as  roughly 


R 


o 


4.251 


(6) 


where  W  is  the  energy  within  the  sphere  in  kilotons  and  pQ  is  the 
ambient  pressure  in  atmospheres .  Rq  is  roughly  the  radius  at  which 
radiative  expansion  gives  way  to  a  shockfront  at  the  periphery  of 
the  fireball.  The  properties  of  the  main  shockwave  are  not  strongly 
dependent  on  the  choice  of  initial  sphere  size  (ref.  4).  Deviating 
from  the  above  value  by  even  a  factor  of  5  does  not  affect  the  main 
shock  seriously,  provided  that  one  does  not  use  the  computed  shock- 
wave  values  until  the  initial  sphere  has  about  doubled  its  radius . 
Some  time  is  needed  for  the  energy  to  be  fed  into  the  shockwave. 


The  weapon  mass  is  usually  neglected  in  these  calculations. 
It  is  usually  fairly  small  compared  to  the  mass  of  air  within  the 
initial  sphere.  Actual  details  of  the  weapon  expansion  can  be 
included  if  the  initial  configurations  and  equations  of  state  for 
the  weapon  materials  are  incorporated. 

The  energy  to  put  into  the  initial  hot-air  sphere  is  Just 


E(erg/gram)  -  k.185  x  101’  w  r — i -  *  E0  (7) 

3 "  V  p 

where  W  is  again  the  energy  in  kilotons  and  p  is  the  air  density 
within  the  fireball  (usually  assumed  to  equal  ambient  density). 

E0  is  the  energy  of  the  uibient  air  before  the  explosion j  it  is  usually 
negligible. 

The  continuous  actual  problem  must  be  divided  into  zones 
for  the  finite  difference  solution.  If  these  zones  are  made  small, 
the  problem  will  run  slowly.  If  the  zones  are  made  larger  to  make 
the  problem  run  faster,  the  results  will  be  less  accurate.  When 
the  zones  are  too  large,  the  shock  tends  to  have  too  high  a  pressure 
and  to  travel  too  fast.  A  useful  compromise  between  computing  time 
and  accuracy  has  been  found  to  be  a  zone  size,  in  the  ambient  air 
ahead  or  the  shock,  of  roughly 


Ax  =  1.0  (W/p0) 


meters. 


(8) 
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3.2  Other  Problems :  This  code  is,  of  course,  not  restricted  to 
nuclear  explosion  calculations.  The  main  adaptations  for  nuclear 
explosion  use  are  the  AKDC  atmosphere,  the  rezoning  routines,  and 
the  summary  routines.  Given  appropriate  equations  of  state,  the 
code  can  be  used  for  a  variety  of  one -dimensional  problems ,  e.g., 
high  explosives  in  vacuum,  air,  water,  or  any  medium  whose  equation 
of  state  is  known. 
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Example:  1  It  E&pioaion  at  See  Level 

The  problem  must  be  divided  into  zones  and  regions  and  initial  values 
must  be  assigned  to  x ,  u,  ©,  and  density  (l/v)  for  each  zo no.  It  is 
convenient  to  divide  this  problem  into  three  unifora  regions  as  shown  in 
the  figure  below:  the  fireball,  a  region  of  fine  senes  for  the  shockwave 
to  move  into,  and  an  outaraost  extensive  region  of  coarse  zones  which  will 
l®  subdivided  into  fine  zones  as  tbs  shockwave  ©ores  ou4 . 

Equation  (6)  gives  a  fireball  radius  of  E0=  4.251  ©store.  Equation 
(8)  gives  one  ester  as  a  suitable  zones  else  in  the  fins  zone  region.  Let 
there  be  ton  fine  zones  per  scares  sons  and  uea  30  fins  zones.  Then  region 
2  is?  30  ©store  thick  and  the  radius  of  its  outer  boundary  is  34.251  ©stars. 
Th©  large  zonae  in  tba  third  region  are  10  ©store  thick.  Using  200  of  them 
would  give  a  problem  boundary  of  2034.251  ssstern,  which  is  far  enough  to 
contain  the  phenomena  of  interest.  The  as  numbers  reed  not  be  so  ©sect; 

35  no  torn  and  2000  cetera  are  close  enough. 

The  initial  internal  energy  in  the  fireball,  from  equation  (7),  is 
1.062x10^  erg/gram.  Tba  internal  energy  of  the  eabient  air,  calculated 
from  e  «  pv/(fr-l),  is  2.068x!09  erg/grara. 

A  sector  of  this  problem  appears  as  follows: 


fixed  outer 
boundary  of 
problem 
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THE  DATA  FOR  THIS  PROBLEM  MIGHT  BE  AS  FOLLOWS- 

CARO  1  WUNDY  CODE  CALCULATION  FOR  1  KT  AT  SEA  LEVEL. 


CARD  2 


CARD  3  NPR08 


CARD  4 


U.S.N.O.L. 


3-1-65 


1  KT  AT  SEA  LEVEL. 


1160 


CARD  5 


CARD  6 


CARD  7 A 


K 

3 

MUR  IN 

2 

MUREX 

2 

IMAX 

3 

I  ARDC 

0 

INTAPE 

0 

INCODS 

0 

NOUIT 

4000 

NPR 

50 

NTAPE 

0 

JCALC 

20 

KOUT2 

4 

KOUT2A 

10 

KOUT2B 

0 

KOUT3 

4 

KOUT3A 

10 

K0UT4 

4000 

KREZ1 

0 

KREZ2 

0 

KREZ3 

10 

KREZ3A 

30 

KREZ3B 

41 

KREZ4 

2 

NC ( 12  )  ALL  Z 

NC( 11) 

1 

NC( 13-24) 

ALL 

NC ( 17 ) 

1 

I  I 

1 

NEQST ( 1 ) 

2 

GAMMA  1(1) 

1.4 

NZONES ( 1 ) 

10 

OUTBDY ( 1  ) 

425 

E I N I  T  (  1  ) 

1.06; 

UINIT(l) 

0. 

D I N  I  T  (  1  ) 

.001, 

C I  NQ  (  1  ) 

2.0 

AINQ( 1  ) 

0.2 

ZONING! 1 ) 

0. 

{ARBITRARY  PROBLEM  NUMBER} 

(SPHERICAL  GEOMETRY) 

(ORIGIN  IS  FIXED) 

(OUTER  BOUNDARY  OF  PROBLEM  IS  FIXED) 

(THREE  REGIONS) 

( ARDC  ATMOSPHERE  NOT  USED) 

(NO  DATA  INPUT  FROM  TAPE  18) 

(NO  EXTRA  DATA  INPUT  FROM  CARDS) 

(STOP  PROBLEM  AT  CYCLE  4000) 

(PRINT  AFTER  EVERY  50  CYCLES) 

(DO  NOT  WRITE  TAPE  18) 

(START  WITH  ONLY  THE  FIRST  20  ZONES  ACTIVE) 

(USE  8  POINTS  TO  GET  EXTRAPOLATED  PRESSURE) 
(STORE  SHOCK  FRONT  DATA  EVERY  10  CYCLES) 

(NO  INTERACTION  BETWEEN  0UT2  AND  OUT3) 
(GATHER  P-T  DATA  AT  4  DISTANCES.  S ( 1-4 ) ) 
(GATHER  P-T  DATA  EVERY  10  CYCLES) 

(PUNCH  CONTINUATION  CARDS  AT  CYCLE  4000) 

( REZ 1  NOT  USED) 

( REZ2  NOT  USED) 

(USE  10  FINE  ZONES  PER  COARSE  ZONE) 

(USE  A  TOTAL  OF  30  FINE  ZONES) 

(OUTER  INTERFACE  OF  LAST  FINE  ZONE) 

(REDUCE  REGION  1  TO  2  ZONES  AT  TIME  TREZO) 


0  EXCEPT- 

(PRINT  OVERPRESSURE 


RATHER  THAN  ABSOLUTE  PR) 


ZERO  EXCEPT- 

( PR  I  NT  INTEGRAL  OF  D(DX)  RATHER  THAN  ZONE  MASS) 

(DATA  FOR  REGION  1  FOLLOWS) 

(USE  REAL  AIR  EQUATION  OF  STATE.  EQS2 ) 

(NOT  USED.  EQS2  COMPUTES  GAMMA  FOR  EACH  ZONE) 
(USE  10  ZONES  IN  THE  FIREBALL) 

,1  (OUTER  RADIUS  OF  REGION  1) 

!E+14  (INTERNAL  ENERGY  OF  MATERIAL  IN  REGION  1) 
(START  WITH  ALL  FIREBALL  ZONES  AT  REST) 

>25  (USE  AMBIENT  SEA  LEVEL  DENSITY) 

(THE  CHOICE  OF  CINQ  AND  A I NQ  IS  AN  ART.  THESE 
VALUES  HAVE  GIVEN  SATISFACTORY  RESULTS) 

(USE  ZONES  OF  UNIFORM  WIDTH) 
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CARD  7B  II 


2  (DATA  FOR  REGION  2  FOLLOWS) 


NEQST ( 2 )  2  (USE  REAL  AIR  EQUATION  OF  STATE*  EQS2 ) 


GAMMA  I { 2 )  1.4  (N 
NZONESt  2 )  30  (U 

OUTBDY ( 2 )  3500. 
EIN I T ( 2 )  2.068E+9 


UINIT ( 2 ) 
DI N I T ( 2 ) 
CINQ ( 2 ) 
AINQ( 2 ) 
20NING( 2 ) 


0. 

.001225 

2.0 

0.2 

0. 


(NOT  USED) 

(USE  30  ZONES  IN  THE  FIRST  AIR  REGION) 

(OUTER  RADIUS  OF  REGION  2  IS  35  METERS) 
+9  (ENERGY  OF  AMBIENT  SEA  LEVEL  AIR) 

(START  AT  REST) 

5  (SEA  LEVEL  AMBIENT  DENSITY) 


CARD  7C  SAME  AS  CARD  7B  EXCEPT- 

II  3  (DATA  FOR  REGION  3  FOLLOWS) 

NZONES ( 3 )  200  (USE  200  ZONES  IN  REGION  3) 

OUTBDY ( 3 )  2.E+5  (OUTER  BOUNDARY  OF  PROBLEM  IS  AT  2000  METERS) 


CARD  8 


COSPHI 

HOBKM 

WKT 

BLANK1 

BLANK2 

ZONSIZ 


(NOT  USED  SINCE  !ARDC=0) 

(CENTER  OF  EXPLOSION  IS  AT  SEA  LEVEL) 
(NOT  USED) 

(NOT  USED) 

(NOT  USED) 

(NOT  USED) 


CARD  9  T 


0.  (START  WITH  TIME=0.) 


TREZO 
DTMI N ( 2  ) 
DTRATE 
STABIL 
UCUT 


(REDUCE  REGION  1  TO  2  ZONES  AT  .01  SECOND) 
(USE  BUILT-IN  VALUE  OF  l.E-10  SECONDS) 

(USE  BUILT-IN  VALUE  OF  1.4) 

(USE  BUILT-IN  VALUE  OF  0.801) 

(USE  BUILT-IN  VALUE  OF  l.E-8  CM/SEC) 


CARD  10  TL I  ST ( 1-6 )  ALL  ZERO  (PRINTING  IS  CONTROLLED  BY  NPR ) 

CARD  11  S(l)  1800.  (COLLECT  P-T  DATA  AT  18.  38*  90*  AND  250  METERS) 


CARD  11  S(l) 
S  ( 2  ) 
S(  3  ) 
S  ( 4 ) 
S(  5) 
S  ( 6 ) 


1800. 

3800. 

9000. 

25000. 

0. 


(NOT  USED) 


0.  (NOT  USED) 


CARD  12  S ( 7-12 )  ALL  ZERO  (NOT  USED) 

CARD  13  BLANK  EXCEPT  FOR  0  IN  COLUMN  5  (NO  MORE  PROBLEMS  FOLLOW) 
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THE  DATA  CAROS  FOR  THIS  PR08LEM  MIGHT  APPEAR  AS  FOLLOWS- 


JDATA 


WUNOY 

COOE 

CALCUIAT 

ION  FOR  1  KT 

AT 

SEA  LEVEL. 

1 

U.S.N. 

O.L. 

3-1-65 

1 

KT  AT 

SEA 

LEVEL 

• 

2 

1160 

3 

2  2 

3 

0 

0 

0 

4000 

50  0 

20 

3 

4 

10 

0  4 

10 

4000 

0 

0 

10 

30  41 

2 

4 

0 

0 

0  0 

0 

0 

0 

0 

0 

0  1 

0 

5 

0 

0 

0  0 

1 

0 

0 

0 

0 

0  0 

0 

6 

1  2 

1.4 

10  4.25100+2 

1.0620+14 

0. 

1.22500-3 

2. 

•  2 

0. 

7A 

2  2 

1.4 

30 

3500. 

2.06800+9 

0. 

1.22500-3 

2. 

•  2 

0. 

7B 

3  2 

1.4 

200  2.00000+5 

2.06786+9 

0* 

1.22500-3 

2. 

.2 

0. 

7C 

0. 

0. 

0. 

0* 

0. 

0. 

8 

0. 

•  01 

0. 

0. 

0* 

0. 

9 

0. 

0. 

0. 

0. 

0. 

0. 

10 

1800. 

3800. 

9000. 

25000. 

0. 

0. 

11 

0 

0. 

0. 

0. 

0. 

0. 

0. 

12 

THE 

RUNNING  TIME 

FOR 

THIS  PROBLEM  IS 

20  MINUTES  ON  AN 

IBM 

7090. 
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Sows  of  the  possible  aodifie&ticns  of  this  problem  are  as  follows: 

1}  The  prograa  will  calculate  its  own  fireball  internal  energy  if 
SINIT(1)»0.  and  WET  is  positive  (hare,  HX?»1.). 

2)  To  calculate  upward  from  the  explosion  rather  than  horizontally,  eat 
IARDC  =2  to  order  the  ARDC  subroutine  to  provide  the  ambient  conditions 
in  regions  2  and  3  and  set  COSPHI'l.  to  indicate  an  upward  ray.  To 
calculate  upward  at  a  45  degree  angle,  set  C0SPHI«,7O7,  ate.  The  ARDC 
subroutine  can  of  course  be  used  for  homogeneous  runs,  if  desired.  If 
IARDC  were  2  in  section  1)  above,  BINIT (2),  EINXT{3),  DINXT(2),  and 
DINIT(3)  would  be  supplied  by  ARDC  and  could  be  zero  in  the  input  cards. 

3)  To  calculate  downward  from  1  ET  at  20  km  altitude  of  burst,  get 
yKT=l. ,  IARDC -2,  C0SFHX=-i. ,  and  H0BKM=20. 

4)  If  tape  18  was  written  during  a  problem  (e.g.,  bv  setting  NTAPfi“500 
to  cause  a  record  to  be  written  every  500  cycles)  the  problem  can  be 
continued  later  from  any  record  of  this  tape.  Use  the  same  data  cards 
as  originally  used  to  run  the  problem  except  for  setting  NQUIT  to  some 
new  cycle  number  to  stop  on  and  INTAPE^  (or  whatever  record  number 

is  desired) . 

5)  To  calculate  one  megaton  rather  than  one  kiloton,  change  the  OUT  BUT 
values  and  WKT .  To  use  the  scaling  criteria  given  here  (Sachs  scaling) 
simply  multiply  the  OUTBDI  values  for  one  kiloton  by  ten.  Actually, 
for  a  homogeneous  calculation  (COSHIIeO.) ,  the  1  KT  results  could  be 
used  directly  for  any  yield  by  applying  the  proper  scaling  factors 

to  all  quantities  in  the  output.  It  is  generally  more  convenient, 
however,  to  re-run  the  problem  for  each  different  yield.  The  results 
for  nonhomogeneity  runs  (COSPKI  non-zero)  noi,  Sachs  scale  and  each 
yield  and  altitude  of  burst  must  be  calculated  separately. 

6)  To  calculate  the  pressure  incident  on  the  ground  from  an  air  burst, 
proceed  as  in  3)  above  but  adjust  the  zoning  so  0UTBDI(3)  equals  the 
distance  to  the  ground.  The  shockwave  will  reflect  from  this  rigid 
boundary  as  from  the  ground,  but  the  problem  will  become  meaningless 
soon  after  reflection  because  the  spherical  symmetry  will  force  the 
shock  to  converge  back  on  the  origin  rather  than  behave  as  a  true 
reflected  shock  apparently  originating  from  a  focus  beneath  the 
ground.  Calculation  of  reflected  shocks  and  Mach  stems  would  require 
a  two-dimensional  code. 
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APPENDIX  8 

LISTINGS  FOR  WUNDY  HYDROCODE 
GLOSSARY  OF  DIMENSIONED  QUANTITIES 


I  =  REGION  INDEX*  RUNS  FROM  1  TO  30  OR  31 
J  *  ZONE  INDEX*  RUNS  FROM  1  TO  300  OR  301 
*  =  NORMALLY  READ  IN  AT  START  OF  PROBLEM 

**  =  NOT  NORMALLY  READ  IN  BUT  CAN  bE  READ  FROM  CARDS  BY  6ENSUB 


ARTIFICIAL  VISCOSITY  TERM* 

TO  BE  PRINTED  AT  START* 

TU  bE  PRINTED  ON  EVERY  OUTPUT  CYCLE. 


*  A I  NO ( I )  COEFFICIENT  IN  LINEAR 

*  AL I  ST ( 1-12 )  ARBITRARY  LETTERING 

*  AL I  ST ( 12-24 )  ARBITRARY  LETTERING 

B ( J )  NOT  USED. 

C(J)  NOT  USED. 

*  CINQ(I)  COEFFICIENT  IN  QUADRATIC  ARTIFICIAL  VISCOSITY  TERM. 

**  D { J 5  DENSITY  (G/CC). 

*  DIN  1 T I  I )  INITIAL  VALUE  OF  D  FOR  ALL  ZONES  IN  REGION  I  (G/CC). 

DTMIN(l)  TIMESTEP  CENTERED  AT  TIME  N-l/2  . 

(2)  TIMESTEP  CENTERED  AT  TIME  N+l/2  . 

(3)  TIMESTEP  CENTERED  AT  TIME  N. 

DTZJ(J)  TIME  STEP  ASSOCIATED  WITH  ZONE  J. 

DUDTIJ)  ACCELERATION  OF  INTERFACE  J. 

DUMMY! J)  AUXILIARY  VARIABLE*  USED  BY  GENR  WHEN  READING  CARDS. 


DWAS( J) 

**  E  C  J ) 

*  EINITU  ) 
EINT(  I  » 
EKIN  (  I  ) 
ENTOT ( I ) 
EWAS ( J ) 

F(  J) 

*  GAMMA I ( I ) 

GAMMA J  {  J ) 
GRAMS! J) 
INTERJ! I ) 
JIN!  J) 

*  NC  < 1-24) 

*  NC ( 1 ) 

*  NC ! 2 ) 

*  NC ( 3  ) 

*  NC { 4 ) 

*  NC ( 5 ) 

*  NC ( 6  ) 

*  NC ( 7  5 


VALUE  OF  D ( J )  AT  PREVIOUS  TIME  STEP. 

INTERNAL  ENERGY  OF  ZONE  J  (ERGS/GRAM). 

INITIAL  TOTAL  INTERNAL  ENERGY  OF  REGION  I 
TOTAL  INTERNAL  ENERGY  OF  REGION  I 
TOTAL  KINETIC  ENERGY  OF  REGION  I 
TOTAL  ENERGY  (INTERNAL  +  KINETIC)  OF  REGION  I 
INTERNAL  ENERGY  E  <  J )  AT  PREVIOUS  TIME  STEP. 

NOT  USED. 

GAMMA  FOR  REGION  I.  IF  GAMMAJ  IS  CALCULATED  FOR  EACH 
ZONE  BY  EQST  *  IT  SUPERCEDES  GAMMA  I • 

GAMMA  FOR  ZONE  J,  (EQUATION  OF  STATE  CONSTANT). 

MASS  OF  ZONE  J  (GRAMS) . 

NUMBER  OF  THE  INTERFACE  FORMING  INNER  BOUNDARY  OF  REG  1. 
AUXILIARY  VARIABLE.  USED  BY  GENR  WHEN  CARDS  ARE  RtAD. 
LOGICAL  FLOW  CONTROL  VARIABLES*  NORMALLY  ZERO. 

=1  LIMIT  THE  RATE  OF  Q  GROWTH  (VERY  RARELY  USED), 

NOT  USED. 

NOT  USED. 

=1  SUPPRESS  ENERGY  CHECK.  FOR  USE  IN  PROBLEMS  WHERE 
ENERGY  IS  NOT  CONSERVED  (E.G.*  MOVING  PISTONS). 

=1  LEAVE  OUT  GRAVITY  TERM  (FOR  TESTING  SIZE  OF  GRAVITY 
EFFECT  ON  NONHOMOGENEITY  PROBLEMS). 

=1  PRINT  NC ( 6 )  TIMES  PER  DECADE  IN  TIME. 

=1  USE  SPECIAL  Q  (FOR  PROBLEMS  WITH  SEVERE  MOTIONS  NEAR 


THE  ORIGIN) , 
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*  NC  ( 8 )  *1  SUPPRESS  TRIGGERING  OF  REZ3  ON  P  CHAN6E. 

*  NC  ( 9  5  *0  PUT  OUT 2  OUTPUT  ON  TAPE  6.  *2.  TAPE  26. 

*  NC(10)  *0  PUT  0UT3  OUTPUT  ON  TAPE  6.  *2.  TAPE  26. 

*  NC(ll)  =1  PRINT  OVERPRESSURE  RATHER  THAN  ABSOLUTE  PS1. 

*  NC  t 12 )  «1  DUMP  THE  MEMORY  CORES  AFTER  PROBLEM  IS  DONE. 

*  NC  ( 1 3  )  =1  CALL  SUBROUTINE  FIXX  EVERY  CYCLE. 

*  NC ( 14 )  *0  PUT  REZONING  COMMENTS  ON  TAPE  6.  *2.  ON  TAPE  26. 

*  NC ( 15 )  =  USE  NC < 15)  STEPS  OF  ONE  CYCLE  IN  0UT3. 

*  NC ( 16 )  =1  PRINT  BEFORE  AND  AFTER  EACH  TIME  REZ3  IS  USED. 

*  NCI  17 )  =1  PRINT  CUMULATIVE  GRAMS/SQ  CM  INSTEAD  OF  PDYN. 

*  NC  ( 18 )  *  DO  NOT  CARRY  MAX  Q  SEARCH  IN  0UT2  INTO  REGION  NCU8). 

*  NC ( 19 )  *1  PRINT  LIST  OF  ALL  TAPE  18  OUTPUT. 

*  NC{ 20-24)  NOT  USED. 

NC( 25-40)  INTERNAL  FLAGS  FOR  LOGICAL  FLOW  CONTROL. 

NC ( 25 )  *1  CALCULATE  ONLY  THE  TIMESTEP  IN  HYDR  (RARELY  USED). 

NC ( 26  )  NUMBER  OF  ACTIVE  LOCATIONS  IN  0UT3. 

NC ( 27 )  NOT  USED. 

NC( 28-40)  NOT  USED. 

*  NEQST(I)  NUMBER  OF  EQUATION  OF  STATE  TO  BE  USED  IN  REGION  I. 

*  NZONES ( I )  NUMBER  OF  ZONES  IN  REGION  I. 

*  OUTBDY(I)  INITIAL  X-COORDINATE  OF  OUTER  INTERFACE  OF  REGION  I. 

P(J)  PRESSURE  BEFORE  Q  HAS  BEEN  ADDED  (DYNE/SQ  CM). 

PDYN ( J )  DYNAMIC  PRESSURE  (PSI). 

PQ { J )  PRESSURE  (MEGABARS)  =P(J)+Q(J) 

PS  I ( J )  PQ ( J )  IN  PSI  UNITS. 

PWAS(J)  VALUE  OF  P(J)  AT  PREVIOUS  TIME  STEP. 

**  Q ( J )  ARTIFICIAL  VISCOSITY  CONTRIBUTION  TO  PRESSURE. 

*  S(L)  SPARE  INPUT  QUANTITIES.  OUT3  USES  SOME  AS  DISTANCES. 

SOUND ( J )  SOUND  SPEED  IN  ZONE  J. 

TEMP ( J )  TEMPERATURE  (KELVIN). 

*  TL I  ST ( 1-8 )  PRINTING  TIMES.  RARELY  USED.  USUALLY*0. 

PRINT  EVERY  TL I  ST ( 1 )  SECONDS  UNTIL  TLIST(2)  IS  REACHED. 
THEN  PRINT  EVERY  TL 1ST  C  3  >  UNTIL  TLIST(4)  IS  REACHED,  AND 
THEN  STOP  THE  CALCULATION.  WRITE  TV  TAPE  EVERY  TLIST(5) 
SEC  UNTIL  TL 1ST ( 2 )  AND  EVERY  TL 1ST ( 6 )  THEREAFTER. 

**  U(J)  INTERFACE  VELOCITY  (CM/SECOND). 

*  UINIT(I)  INITIAL  VELOCITY  OF  ALL  INTERFACES  OF  REGION  I. 

**  X(JJ  COORDINATE  OF  INTERFACE  (CM). 

XAV(J)  COORDINATE  OF  CENTER  OF  ZONE  (CM). 

Y ( J )  NOT  USED. 

*  Z ( J )  Z ( 1-30 )  ARE  READ  IN  BUT  NOT  PRESENTLY  USED. 

*  ZONING(I)  TYPE  OF  ZONING  IN  REGION  I. 

=0.  USE  ZONES  OF  EQUAL  WIDTH  THROUGHOUT  I. 

•POSITIVE*  LET  THE  RATIO  OF  SIZE  OF  OUTERMOST  ZONE  TO 
INNERMOST  ZONE  IN  I  BE  ZONING! I)  AND  FILL  I  WITH  ZONES  OF 
GEOMETRICAL  PROGRESSION  IN  WIDTH. 
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GLOSSARY  OF  NON-D JMENSIONED  QUANTITIES 


*  *  NORMALLY  READ  IN  AT  START  OF  PROBLEM 


AREA 

*  BLANK1 

*  BLANK2 
BLANKS 

*  COSPHI 


*  DTRATE 
E INSUM 
EKSUM 
ESTART 
ETOTAL 
GACCEL 

*  HOBKM 
I 

*  IARDC 

I  DENT 

*  IMAX 

*  INCODS 


*  INTAPE 
J 

*  JCALC 

JCALC1 

JCALC2 

JDT 

JLAST 

JLAST1 

JMAX 

JMIN 

JMH 

JPH 

JUTEST 

*  K 


*  KOUT2 

*  KOUT2A 

*  KOUT2B 


AREA  OF  INTERFACE  J  (SQUARE  CENTIMETERS). 

NOT  USED. 

NOT  USED. 

NOT  USED. 

COSINE  OF  ANGLE  GIVING  DIRECTION  OF  RAY  BEING  CALCULATED. 
MEASURED  FROM  ZENITH.  +1.  FOR  UP.  0.  FOR  HORIZONTAL. 

-1.  FOR  DOWN.  ETC.  USED  FOR  MULTIPLYING  GRAVITY  TERM  IN 
NONHOMOGENEOUS  ATMOSPHERE  CALCS  FOR  GETTING  ZONE  ALTITUDE. 
RATE  AT  WHICH  TIMESTEP  IS  ALLOWED  TO  BUILD  UP  (USU.  1.4). 
TOTAL  INTERNAL  ENERGY  IN  ENTIRE  PROBLEM. 

TOTAL  KINETIC  ENERGY  IN  ENTIRE  PROBLEM. 

TOTAL  ENERGY  IN  ENTIRE  PROBLEM  AT  BEGINNING  (ERGS). 

TOTAL  ENERGY  IN  ENTIRE  PROBLEM  AT  CURRENT  TIME  (ERGS). 

COMPONENT  OF  GRAVITY  ALONG  RAY  BEING  CALCULATED! CM/SECSQ) • 
ALTITUDE  OF  CENTER  OF  PROBLEM  (KILOMETERS). 

REGION  INDEX,  RANGES  FROM  1  TO  IMAX. 

SET  UP  WITH  ARDC  ATMOSPHERE  IN  REGIONS  IARDC  TO  IMAX. 

=0  DO  NOT  USE  ARDC  AT  ALL. 

USED  FOR  IDENTIFYING  SUBROUTINES  WHEN  BEGINNING  FROBLEM. 
NUMBER  OF  REGIONS  IN  PROBLEM  (UP  TO  30). 

=0  NO  EXTRA  DATA  CARDS  TO  BE  READ  IN  GENR. 

=1  READ  TYPE  I  CARDS. 

=2  READ  TYPE  II  CARDS. 

READ  INITIAL  CONDITIONS  FROM  RECORD  INTAPE  OF  TAPE  18. 

ZONE  INDEX.  USED  IN  SEVERAL  PLACES. 

NUMBER  OF  LAST  INTERFACE  CURRENTLY  BEING  CALCULATED  IN  EACH 
SWEEP.  EVERYTHING  BEYOND  JCALC  IS  ASSUMED  TO  BE  INACTIVE. 

* JCALC-1 
= JCALC-2 

INDEX  OF  ZONE  CONTAINING  SMALLEST  TIME  STEP. 

NUMBER  OF  THE  LAST  INTERFACE  IN  THE  PROBLEM.  (UP  TO  301). 
=JLAST-1 

OUTER  INTERFACE  OF  LAST  ZONE  OF  REGION  BEING  CALCULATED. 

OUTER  INTERFACE  OF  FIRST  ZONE  OF  REGION  BEING  CAl CULATED. 

INDEX  OF  CENTER  OF  ZONE  WHOSE  OUTER  INTERFACE  IS  J. 

INDEX  OF  CENTER  OF  ZONE  WHOSE  INNER  INTERFACE  IS  J. 

SOME  INTERFACE  OUTSIDE  THE  INITIALLY  ACTIVE  REGION. 

PROBLEM  SYMMETRY. 

=1  PLANE 
=2  CYLINDRICAL 
*=3  SPHERICAL 

=0  DO  NOT  CALL  OUT2.  (USUALLY  4  WHEN  OUT2  IS  USED.) 

=L  PREPARE  A  LINE  OF  SUMMARY  EVERY  L-TH  CYCLE. 

=-l  READ  IN  TRAJECTORY  DATA  RX.RY, SPEED  IN  OUT3. 

*0  USE  S(L>  LIST  AS  R ( L )  DISTANCES  IN  OUTS. 

=L  (1  TO  8)  LET  OUT 2  LOCATE  FOR  OUT3  THE  R(L)  WHERE  DESIRED 
OVERPRESSURE  LEVELS  S(L)  OCCUR  AND  LET  OUT3  ACCUMULATE 
P-T  DATA  AT  THESE  DISTANCES. 
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*  K0UT3 

*  X0UT3A 

*  K0UT4 

*  KREZl 

*  KREZ2 

*  KREZ3 

*  KREZ3A 

*  KREZ3B 

*  KREZ4 

KUTOFF 

L 

MAS 


*  MUREX 


*  MURIN 


N 

*  NCYCLE 
NDOALL 
NGEOM 

*  NPR 
NPRINT 

*  NPROB 
NQST 

*  NQUIT 


NRECRO 

*  NTAPE 

NWRITE 

*  STABIL 

*  T 

*  TREZO 

*  UCUT 

USQ 

VOL 


NUMBER  OF  DISTANCES  (OR  OVERPRESSURES*  IF  NC(8)  POSITIVE) 

IN  OUT 3*  OUTS  IS  NOT  USED  IF  KQUT3-0. 

STORE  A  LINE  OF  DATA  FOR  EACH  DISTANCE  EVERY  KOUT3A  CYCLES. 
«1  USE  0UT4. 

«1  USE  OUT1. 

»L  USE  REZ2  BEGINNING  WITH  REGION  L. 

■0  OMIT  REZ3 »  *L  USE  L  FINE  ZONES  PER  LARGE  ZONE. 

«  TOTAL  NUMBER  OF  FINE  ZONES. 

J  INDEX  OF  INTERFACE  BOUNDING  OUTERMOST  FINE  ZONE. 

*0  DO  NOT  USE  REZ4.  IF  NON-ZERO*  REDUCE  REGION  1  TO  K.REZ4 
ZONES  WHEN  TIME  TREZO  IS  REACHED. 

«1  TROUBLE.  PRINT  AND  THEN  STOP  THE  PROBLEM. 

GENERAL  SUBSCRIPT.  USED  IN  MANY  PLACES. 

*0  NORMAL  PRINTOUT. 

«1  PRINT  ZOHE  MASSES  INSTEAD  OF  DYNPSI  (USED  INITIALLY 
AND  AFTER  REZONING). 

TYPE  OF  OUTER  BOUNDARY  {INTERFACE  JLAST )  ON  PROBLEM. 

«1  FREE  SURFACE. 

■  2  RIGID  WALL. 

*3  BOUNDARY  MOTION  TO  BE  SPECIFIED  BY  SUBROUTINE  BDY2. 

TYPE  OF  INNER  BOUNDARY  (INTERFACE  1)  ON  PROBLEM. 

*1  FREE  SURFACE. 

*2  RIGID  WALL. 

*3  BOUNDARY  MOTION  TO  BE  SPECIFIED  BY  SUBROUTINE  BDYl. 

NOT  USED. 

CYCLE  NUMBER.  EACH  CYCLE  IS  ONE  COMPLETE  SWEEP  THRU  MESH. 
*1  CALCULATE  ALL  ZONES  WHETHER  ACTIVE  OR  NOT. 

SPECIFIES  WHAT  IS  WANTED  FROM  SUBROUTINE  GEOM  (Q.V.). 

«1  CALCULATE  AREA,  =2  CALCULATE  VOLUME. 

PRINT  OUTPUT  EVERY  NPR  CYCLES  (USUALLY  ABOUT  25). 

*1  FORCES  PRINTING  AFTER  CURRENT  CYCLE. 

ARBITRARY  PROBLEM  NUM8ER. 

“NEQST ( I ) . 

TERMINATE  THE  PROBLEM  WHEN  NCYCLE  REACHES  NQUIT. 

(WUNDY  WILL  DO  ABOUT  2000  CYCLES  FOR  40  ACTIVE  ZONES  IN 
ABOUT  10  MINUTES  ON  THE  7090) 

NUMBER  OF  RECORD  ON  TAPE  18.  (1-FIRST  RECORD*  ETC.). 

WRITE  TAPE  18  EVERY  NTAPE  CYCLES  (USUALLY  500). 

»0  DO  NOT  USE  TAPE  18. 

FORCES  WRITING  OF  TAPE  18  AFTER  CURRENT  CYCLE. 

STABILITY  CONSTANT  (USUALLY  ABOUT  0.8). 

CURRENT  TIME  (SECONDS). 

CALL  REZ4  WHEN  TIME  TREZO  IS  REACHED. 

SET  ALL  INTERFACE  VELOCITIES  LESS  THAN  UCUT  EQUAL  TO  0. 

USE  UCUT *l.E-2  IF  UCUT  =  0.  IN  THE  INPUT. 

SQUARE  OF  INTERFACE  VELOCITY. 

VOLUME  OF  ZONE  BOUNDED  BY  INTERFACES  J  AND  J-l  . 
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*  WKT  FIREBALL  ENERGY*  K1LOTONS.  IF  THIS  IS  NON-ZERO*  PROGRAM 

SUPPLIES  VALUtS  FOR  ANY  OF  THE  FOLLOWING  IF  THEY  ARE 
ZERO  IN  THE  INPUT  — 

N  ZONES  (  1-3 )  » OUT  BOY  (  1-3}  *TREZO»DINIT  (1)  *EINITU1 

*  ZONSIZ  AMBIENT  ZONE  SIZE  (METERS)  FOR  1  K.T  AT  SEA  LEVEL*  USED 

BY  PROGRAM  WHEN  IT  SETS  UP  ITS  OWN  ZONING  BY  SACHS 
SCALING.  (PROGRAM  USES  .5  IF  ZONSIZ=0.). 


LIST  OF  SUBROUTINES 


MAIN 
GENR 
HYDR 
OUT  1 
OUT2 
OUT  3 
OUT  A 
REZ1 
REZ2 
REZ3 
REZ4 
ARDC 
GEOM 
EQST 
EQS1 
EOS2 

**  EQS3-12 
**  BDY1 
**  BDY2 
**  F 1 X  X 


HANDLES  MAIN  LOGICAL  FLOW. 

READS  INPUT  AND  INITIALIZES  THE  PROBLEM, 

COMPUTES  HYDRODYNAMIC  MOTIONS  AND  NEXT  TIMESTEP. 

PRINTS  NORMAL  OUTPUT  ON  TAPE  6  AND  WRITES  TAPE  18. 
ACCUMULATES  DATA  ON  MAIN  SHOCK.  FRONT  POSITION*  ETC. 
ACCUMULATES  PKESSURt.  ETC  VS  TIMt  DATA  AT  FIXED  POSITIONS. 
PUNCHES  CARDS  FROM  WHICH  THE  PROBLEM  MAY  BE  CONTINUED. 
REMOVES  EXCESSIVELY  COMPRESSED  ZONES. 

EXPANDS  RANGE  OF  PROBLEM  oY  FACTOR  OF  TWO. 

KEEPS  FINE  ZONES  IN  SHOCK. 

REZONES  FIRST  REGION  AT  TIME  TREZO. 

CONTAINS  ARDC  MODEL  ATMOSPHERE. 

DESCRIBES  PRObLEM  GEOMETRY. 

CONTROLS  ACCESS  TO  EQUATIONS  OF  STATE. 

EQUATION  OF  STATE  FOR  GAMMA  LAW  GAS*  E=PV/(G-1). 
APPROXIMATE  EQUATION  OF  STATE  FOR  REAL  AIR. 

RESERVED  FOR  OTHER  EQUATIONS  OF  STATE. 

SPECIFIES  MOTION  OF  INNER  INTERFACE. 

SPECIFIES  MOTION  OF  OUTER  INTERFACE. 

USED  FOR  MAKING  SPECIAL  CHANGES  DURING  A  RUN. 


**  NOT  INCLUDED  IN  THESE  LISTINGS. 


UNUSED  SUBROUTINES  MAY  BE  REPLACED  BY  DUMMY  DECKS*  E.G.* 

$  I BFTC  BDY2  LIST  *REF 
SUBROUTINE  BDY2 
NULL=0 
RETURN 
END 
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STANDARD  COMMON  STATEMENT 


COMMON  A !NQ( 30)  »ALI5T<24)  »B(300)  .C5300) 

1  »CINO(  30  )  »D(  300 )  »DINI T  <  30 )  .DTMIN(3)  »DTZJ(.300) 

2  *DUDT  ( 301 )  » DUMMY  { 301 )  .DWASOOO)  »E(300)  .EINITJ30) 

3  tEINT ( 30 )  »  E  K I N  < 30)  .ENTOT(30)  »EWAS(300)  .FJ300) 

4. GAMMAK30)  »GAMMA  J  ( 300  )  *GRAMS  ( 300  )  »INTERJ{31)  .JINC301) 

5. NC(40>  *NEQST  { 30 )  .NZONES(30)  .OUTBDYOO)  *P(300) 

6. PDYNC300)  »PQ( 300 )  .PSIC300J  »PWAS(300)  »Q(300) 

7.S ( 24 )  .SOUND(300)  .TEMP<300)  .TLIST(8)  *U ( 301 ) 

8»UI NI T ( 30 )  .X ( 301 )  .XAV<300>  • Y ( 301 )  »Z( 301 ) *ZONING( 30) 

COMMON  AREA  .BLANK1 .BLANK2 ♦ BLANK3 .COSPHI .DTRATE ♦ EINSUM. EKSUM 

1  tESTART  tETOTAL  .GACCEL  .HOBICM  ,  I  , IARDC  *  IDENT  »IMAX  *  I  NCODS 

2  » INTAPE* J  * JCALC  # JCALC1* JCALC2  * JDT  » JLAST  .JLAST1.JMAX 

3#  JMI N  *  JMH  ♦  JPH  *  JUTEST » X  ♦XOUT2  .KOUT2A .KOUT2E .KOUT3 

4. K0UT3A.K0UT4  .KREZ1  ,KREZ2  .KREZ3  .KREZ3A.KREZ3B.KREZ3C .KREZ4 

5 . KUTOFF »L  .MAS  .MUREX  .MURIN  »N  .NCYCLE.NDOALL .NGEOM 

6. NPR  .NPRINT.NPROB  .NOST  . NQU I T  . NRECRD .NTAPE  .NTVREC .NWRI TE 

7.STABIL.T  .TREZO  .UCUT  »USQ  .VOL  *WKT  .ZONSIZ 


nnn  n  n  n  no  n  *  n  n  n  n  v * 
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1BFTC  MAIN  LIST *REF*DO 

MAIN  PROGRAM  1-D  HYDROCODE  MAIN 

NOL  ONE-DIMENSIONAL  HYDRODYNAMIC  CODE. 

THIS  PROGRAM  IS  SET  UP  FOR  CM-GRAM-SECOND  UNITS. 

****  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 

1  IDENT  *0 

CALL  GENS 

BEGIN  IDENTIFICATION  SECTION. 

WRITE (6, 5)  * 

a  FORMAT ( 4QH1MAIN  PROGRAM  VERSION  3.  3-  1-65  ) 

IDENT=1 

CALL  GENR 
CALL  HYDR 
CALL  OUT  1 

I F { KOUT2 .GT .0 )  CALL  OUT2 
I F ( KOUT3 • GT  .0 )  CALL  OUT3 
I F ( K0UT4. GT .0 )  CALL  OUT4 
IF(KREZI.GT.O)  CALL  REZ1 
IF (KREZ2.GT.0I  CALL  REZ2 
IFIKREZ3.GT .0)  CALL  REZ3 
IF ( KREZ4.GT .0 )  CALL  REZ4 
I F ( NC ( 13 ) .GT. 0 )  CALL  FIXX 
CALL  GEOM 
20  DO  22  I *1 » I  MAX 
NQST=NEQST ( I  ) 

CALL  EOST 
22  CONTINUE 
IDENT=0 

END  IDENTIFICATION  SECTION. 

PRINT  INITIAL  CONDITIONS. 

NPRINT=1 
MAS=1 

CALL  OUT  1 

BEGIN  MAIN  CALCULATION  LOOP. 

CHECK  IF  PROBLEM  IS  COMPLETED  YET. 

25  I F ( T »GE. TL I  ST ( 4 ) )  GO  TO  50 

26  IF(NQUIT)  50*28.27 

27  I F ( NCYCLE. GE.NQU I T )  GO  TO  50 

C  STOP  PROBLEM  IF  ZERO  TIMESTEP  OCCURS. 

28  IF(DTMIN(2).GT.O.  )  GO  TO  31 

29  WR I TE ( 6 *30 ) 

30  FORMAT ( 18 HO TIMESTEP  IS  BAD.  ) 

GO  TO  50 

31  T=T+DTMIN12) 

NCYCLE*NCYCLE+1 

CALL  HYDR 


l 
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CALL  OUT 1  MAIN 

RESET  CALCULATE-EVERY-ZONE-REGARDLESS  INDICATOR. 

REMOVE  THE  FOLLOWING  STATEMENT  IF  ALL  ZONES  ARE  TO  BE  CALCULATED 
THROUGHOUT  THE  PROBLEM. 

ND0ALL*0 

IF(KOUTZ.GT.O)  CALL  OUT2 
IF ( KOUT3.GT .0 }  CALL  OUT3 
IF ( K0UT4.GT »0 )  CALL  OUT4 
I F ( KREZ1 sGT.O )  CALL  REZ1 
I F ( KREZ2.GT .0 )  CALL  REZ2 
I F  ( ICREZ3.GT .0  )  CALL  REZ3 
I F ( KREZ4.GT .0 )  CALL  REZ4 
I F ( NC ( 13 ) *GT»0 )  CALL  FIXX 
GO  TO  25 

END  MAIN  CALCULATION  LOOP. 


50  WRITE(6»51)  * 

51  FORMAT ( 40H0  NORMAL  END  REACHED  IN  MAIN  ROUTINE  ) 

NPRINT=1 

NWR I T£*l 

CALL  OUT  1 


C 

C 


C 


SETTING  NQUIT  NEGATIVE  CAUSES  FINAL  PRINTOUTS, 

NQUITs-1 

IF ( KOUT2.GT .0 )  CALL  OUT2 
IFIKOUT3.GT.O)  CALL  OUT3 
IF (KOUT4.GT .0 )  CALL  0UT4 
IF(NC(13).GT.O)  CALL  FIXX 

WRITE  END-OF-FILE  ON  TAPE  26  (IF  USED), 
IF(NC(9).EQ.2.CR,NC(10),EQ,2)  END  FILE  26 
PRINT  TERMINATION  CONDITIONS, 

57  WR  I  TE  ( 6 » 58  )  * 

58  FORMAT ( 3 3H1 PROBLEM  CONDITIONS  AT  END  OF  RUN) 

WR I TE ( 6 » 59 )  I  MAX .KREZ3B , JCALC  »NRECRD  * 

59  FORMAT ( 6H0 I MAX= 1 2  » 10H  KREZ3B  =  I 3 . 10H  JCALC*I3,10H  NRECRD* I  3 ) 

WR I TE ( 6  *60 )  ( N ZONES (L) » L  =  1 » I MAX)  * 

WR I TE ( 6  *  63 )  NCYCLE.T  * 

63  FORMAT (8H0NCYCLE=I5»8H  T IME= 1PE 1 1 . 3 > 

60  FORMAT ( 8H0NZ0NES= 15 14 ) 

DO  ANOTHER  PROBLEM  IF  NEXT=1 

61  READ (5*62)  NEXT  * 

62  FORMAT ( 15) 

99  I F ( NC ( 1 2 ) ,GT • 0 )  CALL  PDUMP ( A I NQ( 1 ) , ZONS IZ » 1 ♦ AI NO ( 1 ) »ZONS IZ »2 ) 

I F ( NEXT • EQ» 1 )  GO  TO  1 
STOP 
END 
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t 


$ IBFTC  GENR  LIST *R£F 

SUBROUTINE  GENR 

INPUT  GENERATOR  GENR 


*****  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 
C 

1  FORMAT ( 12A6  ) 

2  FORMAT (1215) 

3  FORMAT ( I2.I3.F6.3.I4.4E10.5.3F5.2) 

4  FORMAT I 6E1Q.4 ) 


5  FORMAT { 53H1 

PRINTOUT 

OF  INITIAL  CONSTANTS 

) 

6  FORMAT  I 1H0 ♦ 

12A6 ) 

8  FORMAT (96H0 

NPROB 

K  MURIN  MUREX 

I  MAX 

IARDC 

INTAPE 

1  INCODS 

NQUIT 

NPR  NTAPE  JCALC  ) 

9  FORMAT (1218) 

10  FORMAT ( 96H0 

K0UT2  K0UT2A  K0UT2B  K0UT3 

KOUT3A 

K0UT4 

•CREZ1 

1  KREZ2 

KREZ3  KREZ3A  KREZ3B  KREZ4 

> 

11  FORMAT (96H0 

NC<1) 

NC I  2 )  NC I  3 )  NCI4) 

NCI  5) 

NCI  6) 

NCI  7) 

1  NC  ( 8  ) 

NC ( 9 )  NC ( 10  5  NCI  11 )  NC 1 12 ) ) 

12  FORMAT ( 96H0 

NC  1 13 )  NC 1 14 )  NC 1 15 )  NCI16) 

NCI  17 ) 

NC 1 18 ) 

NCI  19) 

1  NC ( 20 )  NCI  21 )  NC( 

22)  NC I  23 )  NC I  24 ) ) 

13  FORMAT ( 105H0I  NEQST 

GAMMA  NZONES  OUTBDY 

EINIT 

1UINIT 

DINIT 

CINQ  AINQ  ZONING 

) 

14  FORMAT! 12.17. F7. 3. 18. 

1P4E14.5.0P3F6.2) 

15  FORMAT (96H0 

COSPHI 

HOBKM 

WKT 

BLA 

1NK1 

BLANK2 

ZONSIZ 

) 

16  FORMAT ( 96H0 

T 

TREZO 

DTMINI 

2) 

DTR 

1ATE 

STABIL 

UCUT 

) 

17  FORMAT ( 1P6E16.5 ) 

18  FORMAT ( 10H0S ( L )  LIST) 

19  FORMAT ( 15HOPRI NT  I NG  TIMES) 

20  FORMAT(5( I4.E10.3) ) 

C 

I F ( IDENT.NE.l)  GO  TO  30 

28  WR I TE ( 6  *29  )  * 

29  FORMAT (40H0SUBROUT I NE  GENR  VERSION  3,  3-  1-65  ) 

RETURN 

C 


30  READ ( 5  » 1 ) 
READ( 5 » 1 ) 
READ! 5*2) 

1 

READ! 5*2 ) 

1 

READ( 5*2) 
READ( 5*2  ) 
READ( 5*3) 

1 

READ( 5*4) 
READ (5*4) 
READ( 5*4 ) 


(ALIST(L) *L=1.12)  * 
l AL 1ST  I L ) »L=13»24)  * 
NPROB  .K  *MURIN  .MUREX  *IMAX  *IARDC  »  * 
INTAPE. INCODS. NQUIT  ,NPR  ,NTAPE  ,JCALC  * 
K0UT2  .KOUT2A.KOUT2B.KOUT3  .KOUT 3A .K0UT4  *  * 
KREZ1  .KREZ2  .KREZ3  .KREZ3A ,KREZ 3B .KREZ4  * 
(  NC  (  L  )  *L=  1 . 12  )  * 
( NC ( L ) »L= 13  *24 )  * 


(  I  I  .NEQST  (  I  >  .GAMMA  I  (  I  )  *NZONESU  )  .OUTBDYI  I)»EINIT(I).  * 
U IN  IT (I), DINIT (I) ,CINQ{ I ) ,AINQ( I > .ZONINGII) »I*1»IMAX)* 


COSPHI ,H0BKM,WKT  .  BLANK  1 .3LANK2 * ZONS I Z  * 

T.TREZ0,DTMIN(2) .DTRATE .STAB  I L .UCUT  * 


(TLIST(I).I=1,6) 


* 
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C 


C 


READ! 5,4)  (S(L)»L*1,125  * 

KUMBAK»1  GENR 

WRI TE ( 6 » 5 )  * 

31  WR1 TE( 6*6)  (ALIST(L) ,L»1,12)  • 

WRITE (6*6)  <ALIST<L) ,L*13,24)  « 

WR I TE ! 6,8 )  ft 

WR I TE ( 6  *9 )  NPROB  ,K  ,MURIN  .MUREX  »IMAX  .IARDC  »  « 

1  INTAPE, INCODS, NQUIT  ,NPR  »NTAPE  .JCALC  « 

WRITE(6,10)  * 

WR  I  TE ( 6 *9 )  KOUT2  .KOUT2A .KOUT2B  ,KOUT3  ,K0UT3A»KQUT4  ,  * 

1  KREZ1  »KREZ2  »KREZ3  .KREZ3A.KREZ3B.KREZ4  * 

WRITE! 6,11)  * 

WRITE(6,9)  ( NC ( L ) , L* 1 , 12 )  * 

WR I TE ( 6, 12  )  » 

WRI  TE(  6,9)  (NCtL) ,L=13.24)  * 

WRI TE(6,13 )  * 

WR  I  TE  (  6, 14  )  (I  ,NEQST(I  )  .GAMMAI  {  I)  ,NZONESm  .OUTBDYII )  ,EINIT<  I  )  ,  * 

1  UINIT1 I ),DINIT( I ) ,CINQ(I ) ,AINQ( I ) .ZONING! I ) ,I«1,IMAX!« 

WRITEI6.15)  * 

PRINT  ZERO  INSTEAD  OF  HOBKM  AND  WKT. 

ZER0*-0. 

WRITEI6.17)  COSPHI .ZERO  .ZERO  .BLANK1 .BLANK2 .ZONSIZ  * 

WRITE !6,16)  * 

WRI TE ( 6 » 17 )  T,TREZO.DTMIN<2) .DTRATE.STABIL.UCUT  * 

WRITEI6.19)  * 

WR I TE  <  6 , 17 )  ! TL 1ST ! L  J »L  =  1 ,6 )  * 

WRITE16.10)  # 

WRITE(6,17)  ( S ! L ) , L= 1 , 12 )  « 

60  TO  (32,83) .KUMBAK 


MAKE  VARIOUS  ADJUSTMENTS  OF  THE  INPUT. 
VELOCITIES  LESS  THAN  UCUT  WILL  BE  SET  TO  ZERO. 
32  IFtUCUT.LE.O. )  UCUT=.01 

IF (DTMIN! 2 ) .LE.O. )  DTMI N ( 2 ) «1 . E-10 
36  DTMlN(3)»DTMIN(2)/2. 

IF ( DTRATE.LE.O.  )  DTRATE=1 .4 
IF (STABIL.LE.O. )  STABIL«.801 
IF! IARDC.LE.O)  IARDC=31 


THE  SECTION  FROM  475  THRU  501  IS  USED  ONLY  FOR  AUTOMATICALLY 
SETTING  UP  UNIFORM  HOT-AIR-SPHERE  EXPLOSION  CALCULATIONS. 

IT  IS  BYPASSED  IF  WKT*0. 

WKT  IS  YIELD  IN  KILOTONS. 

ZONSIZ  IS  ZONE  SIZE  FOR  1KT/SL  AND  IS  ABOUT  .5  TO  1. 

THE  FOLLOWING  SECTION  (475  THRU  501)  REPLACES  ZEROES  IN  NZONES, 
OUTBDY,  EINIT,  AND  DINIT(l)  BY  SACHS  SCALED  VALUES. 

475  IF (WKT. LE.O. )  GO  TO  50 
479  H0BCM*H0BKM*1 , E+5 

CALL  ARDC ( HOBCM.PRAT IO.TRATIO.DRATIO) 
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THIRD*l./3. 

RSCALE* ( 1 • /PRAT  10) **THl RD 
WSCAL£*WKT**TH I RD 

CAMB«PRAT10*1.01325E+6/(DRATIO*1.225E-3*.4) 

I F ( ZQNS I Z  »  L  E  •  0  •  )  Z0NSIZ*.5 
481  Z51ZE*ZONSIZ*100,*RSCALE*WSCALE 

IF (0UT8DY ( 1 )  .LE.O* J  OUTBDy (1 ) *425. 1*RSCALE*WSCAL£ 
iF(NZONESd).LE.O)  NZONES{l)*iO 
l F ( NZONES ( 2  )  «LE .0 )  NZON£S<2)=30 
ZN2*NZ0NES ( 2 ) 

I F (OUTBDY ( 2 ) .LE.O. >  OUTBDY ( 2 ) *ZS lZE*ZN2f OUTBDY { 1 ) 

488  NZ3LIM=299-NZONES(l )-NZ0NES(2 ) 

I F ( NZONES ( 3 ) *GT .0 )  GO  TO  4890 

489  NZONES(3)=NZ3LIM 
GO  TO  490 

4890  1F(OUTBDY(3).GT.O. )  GO  TO  496 

490  ZN3*NZ0NES ( 3 ) 

ZNELL*KREZ3 

DX3=ZSIZE*ZNELL 

I F ( OUTBDY ( 3 ) .GT .0. )  GO  TO  492 

491  OUTBDY  <  3  I =DX3*ZN3+OUTBDY ( 2 ) 

GO  TO  496 

492  SI ZE3=0UTBDY ( 3  J -OUTBDY ( 2 ) 

NZ3=SIZE3/DX3 

I F (NZ3.LE.NZ3L IM)  GO  TO  496 

494  NZONES ( 3  J  =NZ3 

WR I TE ( 6 ,495 )  NZ3.NZ3LIM 

495  FORMAT ( 33H  REGION  3  HAD  TO  BE  ADJUSTED  FROM, 14, 3H  T0,6H  ZONES) 

496  IF( TREZO.LE.O. )  TREZO= .01*RSCALE*WSCALE 

IF (DINIT (1) .LE.O. )  DINIT(1)=DRAT 10*. 001225 
501  I F  C  E IN  1 T  C 1 S.GT.O. )  GO  TO  50 

EINITU )=WXT*4.185E19/(4.1887902*0UTBDY(1)**3«DINIT(1) )  +EAMB 
END  SPECIAL  INPUT  CALC  SECTION. 

LEAVE  OUT  GRAVITY  IF  NC(5)  IS  POSITIVE. 

50  GACCEL=0. 

IF (NCI  5) •LE.O)  GACCEL=COSPH 1*980. 616 

53  DO  54  L=25.40 

54  NC ( L ) *0 
NCYCLE=0 

58  NRECRD=0 
JDT  =  0 
XUTOFF=0 

C  FORCE  CALC  OF  ALL  ZONES  ON  FIRST  PASS  IN  HYDR. 

NDOALL* 1 
MAS=  1 

IF( TLIST( 1 J.EQ.O. )  TL I  ST ( 1 ) *1 • E20 
I F ( TL I  ST (4) .GT.O. )  GO  TO  61 

59  DO  60  L=2»6 


GENR 
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60  TLIST(L)=1.E20  6ENR 

61  TLIST(7)=T+TLIST(1) 

CALCULATE  REGION  INTERFACE  NUMBERS. 

62  INTERJ ( 1 ) = 1 
DO  63  1=1*1  MAX 

INTER J{ I +1 )= INTERJ ( I ) +NZONES ( I ) 

63  CONTINUE 

INITIALIZE  VARIOUS  QUANTITIES. 

JLAST=INTERJ( IMAX+1 ) 

JLAST1=JLAST-1 
IF( JCALC.EQ.OI  JCALC= JLAST 
JCALC1=JCALC--1 
xm  =0. 

U( 1 )=0. 

Q( 1 )*0. 

DTZ J ( l ) =0 . 

CALCULATE  INITIAL  D  AND  t  FOR  EACH  ZONE 
DO  76  1=1*1  MAX 
JMI N= I NTER J ( I ) +1 
JMAX= INTERJ! I +1 ) 

64  IFIZONINGI I ) )  67,67,65 

65  X  I N  =  OUT0DY ( 1-1 ) 

XOUT  =  OUTBDY ( I ) 

LL  =  NZONES<  n-1 
ZLL=LL 

ZONFAC  =  ZON I NG I  I )**( l./ZLL) 

SUM= 1 , 

TERM= 1 , 

DO  66  L= 1 , LL 
TERM=TERM*ZONFAC 

66  SUM=  SUM+TERM 
ZS I ZE= ( XOUT-X I N ) /SUM 
DELFAC= 1 . 

GO  TO  70 

67  DNZONE=NZONES ( I ) 

I F ( I.LE.l  )  GO  TO  69 

68  DELX= (OUTBDY ( I )-OUTBDY( 1-1 J I/DNZONE 
GO  TO  70 

69  DELX= ( OUTBDY ( 1 ) -X ( 1 1 ) /DNZONE 

C  BEGIN  LOOP  FOR  ZONES  WITHIN  REGION  I. 

70  DO  75  J  =  JM I N » JMAX 
JMH= J-l 

U( J)=UINITI I ) 

IFIZONINGI II.LE.O. >  GO  TO  72 

71  X( J)=X( J-l )+DELFAC*ZSIZE 
DELFAC=DELFAC*ZONFAC 


w 
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GO  TO  720  GENR 

72  X( J)=X( J-l J+DELX 
720  0(J)=0. 

PJ J)-0. 

DTZJ( J)=0. 

I F ( I .LT. IARDC)  GO  TO  74 

C  USE  AROC  ATMOSPHERE  BEGINNING  WITH  REGION  IAROC 

73  ALTCM=COSPHI#( XlJl+XI J-l ))*.5  +H0BKM*l.E+5 

CALL  ARDC( ALTCM ,PRAT 1 0 , TRAT 10, DRAT  10) 

D ( JMH ) =DRAT 1 0*  «00 1 225 
E( JMH) =PRATI0*1 ,01325E6/D( JMH)/. 4 
GO  TO  75 
C 

74  D( JMH)=DINIT( I  ) 

E ( JMH)=E!NIT( I  ! 

75  CONTINUE 

76  CONTINUE 

BEGIN  AUXILIARY  INPUT  SECTION 

77  I F ( INTAPE.LF.O)  GO  TO  146 

READ  BINARY  TAPE  18  IF  J.NTAPE  IS  POSITIVE. 

START  PROBLEM  FROM  RECORD  NUMBER  INTAPfc  OF  TAPE  18. 

78  NSK I PA= INTAPE-1 
I F ( NSK I  PA . LE »0  >  GO  TO  81 

79  REWIND  18 
DO  80  JJ= 1 ,NSK I  PA 


READ  (18)  * 

80  READ  (18)  * 

81  READ  (18)  T,(DTMIN(L) * L  =  1  * 3 ) .NCYCLE ,NRECRD, JLAST . JCALC , JCALC1 .* 

1  KREZ3B » I  MAX  »NPROB  » 

JC ALC 1  =  JCALC-1 
JLAST1=JLAST-1 

READ  (18)  (NZONES( I ) ♦INTERJl 1  +  1 ) »I=1»IMAX) ,X( JLAST) »U(JLAST)  »  * 

1(X(L)»U(L)»D(L)»E(L)»Q(L),PQ(L) .TEMP ( L ) .DTZ J ( L ) ,L«1,JLAST1 )  * 

N1 8=  1 8 

144  WRITEI6.145)  INTAPE. N18  * 


145  FORMAT( 39H0INITIAL  CONDITIONS  READ  IN  FROM  RECORD, 13. 8H  OF  TAPE, 
113) 

READ  TYPE  1  CARDS  IF  INCODS  =1. 

146  I  F ( INCODS. NE.l )  GO  TO  199 

NAMELIST  /L I  ST  1 /  B ,C ,D , DTM IN , DTZ J ,DUDT .DUMMY ,DHAS , E »EWAS »F , 

1  INTERJ*JIN,NC»NEOST,NZONES,P»PQ»PWAS.O*TEMP,U,X,Y,Z. 

2  BLANK.3, 1  MAX, INCODS, JCALC  ,  JLAST .JUTEST .KREZ3B , 

3  NCYCLE,NPRINT,NPROB,NRECRD,NWRITE,T 
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147  READ ( 5 * 1 1  ST] ) 


READ  TYPE  2  CARDS  IF  !NC0DS*2. 

199  I F ( INCODS. NE. 2 )  GO  TO  810 

200  READ ( 3*201 )  L*X(L}»U(L)»0(L)»E(L)»Q(L) 

^01  FORMAT ( 14, 4E14. 7*1612.5) 

IFIL.LT. 300)  GO  TO  200 
202  WR I TE ( 6  *  203 ) 

*■03  FORMA T(  3 8H0 INITIAL  CONDITIONS  READ  IN  FROM  CARDS) 

END  AUXILIARY  INPUT  SECTION 
C 

810  KUM3AK=2 
WR I TE ( 6  » 82  ) 

82  FORMAT  (  29H1INPUT  DATA  AFTER  ADJUSTMENTS) 

GO  TO  31 

8  3  DO  84  L  =  1 » JLAST  1 
DWAS(L)=D(L) 

EWAS(L)=E!L) 

84  CONTINUE 
C 

00  87  1=1*1  MAX 
1=  I 

NOST=NEOST ( I ) 

JM!N=INTERJ(I)+1 
JMAX= I NTER J ( I + 1 ) 

DO  86  J  =  JM I N  *  JMAX 
J  =  J 
JPH=  J 
JMH=  J- 1 

85  NGE0M=2 

CALL  GEOM 

GRAMS! JMH) =D( JMH )» VOL 

NOST  NEOS^f5  p  (  J  *  ♦£  (  J  )  tTEMP  (  J  )  *GAMMA  J  (  J  >  *  ANT  SOUND(J). 
CALL  EOST 

PO ( JMH ) =P ( JMH ) +0 ( JMH ) 

86  CONTINUE 

87  CONTINUE 

88  E I NSUM=0 • 

EKSUM=0. 

ESTART=0. 

00  90  I =1 , IMAX 
t'INT!  I  )  =0. 

EKIN! I )=0. 

JM I N= I NTER J ( I ) 

JMAX  INTERJ! 1+1 )-l 
DO  89  J=JMIN,JM„X 
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JPH=J 

USO= ( U ( J ) **?  +U ( J+ 1 ) **2 )  *  •  5 
HNTI  I  )  *E  I  NT  I  '  )  +  F  (  JPH  )  *GR  AMS  <  JPH  ) 
EK 1 N ( I  1 =EK I N ( I )  +. 5*USO*GRAMS ( JPH) 

89  CONTINUE 

E  I  NSUM=£ I NSUM+t I  NT (  I  / 

EKS'JM  =  EKSUM+EK  !N(  !  \ 

90  CONTINUE 

ES  TART  =  E I N5UM+EKSUM 
ETOTAL=ESTART 
JCALC1=JCALC-1 
JLAST1=JLAST~1 
JUT  EST  =  JCALC 
RETURN 
END 


GENR 
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$  I BFTC  HYDR  LIST. REF 
SUBROUTINE  HYOR 
C 

#**##  INSERT  STANDARD  COMMON  CARDS  BEFORt  COMPILING  #**#* 

C 

IF( IDENT.NE.l)  GO  TO  203 

201  WR I TE  <  6 .202 ) 

202  FORMAT ( 40H0SUBROUT l N£  HYOR  VERSION  3.  3-  1-65  ) 

KU  =  0 

RETURN 

C 

C  CALCULATE  ONLY  THE  TIMESTEP  IF  NC(25)«1  INOT  NOW  USED). 

203  I F ( NC ( 25 ) *GT.O !  GO  TO  260 
C 

C  CALCULATE  MOTION  OF  FIRST  INTERFACE 

205  J=1 
JPH=  1 

MUR  I N=MUR I N 

GO  TO  (206.207.208) .MURIN 

206  NGEOM= 1 
CALL  GEOM 

DUDT ( 1 ) =-PQ ( JPH ) /GRAMS ( JPH } / 2 •  *  AREA 
GO  TO  210 

207  DUDT ( 1 ) =0. 

GO  TO  210 

208  CALL  BDY1 

210  U< 1 )=U< 1 )+DTMIN( 3)*DUDT( 1) 

X( 1 ) =X( 1 )+DTMIN(2 )  *U  ( 1 ) 

BEGIN  MAIN  HYDRODYNAMIC  LOOP 

211  JDONE=0 
DO  250  I  - 1 . 1 MAX 
1  =  1 

JM I N= I NTERJ ( I ) +1 
JMAX*INTERJ( 1+1 ) 

DO  240  J= JMIN , JMAX 
J*J 
JPH  =  J 
JMH= J- 1 

C  CALCULATE  LAST  INTERFACE  SEPARATELY 

I F ( J.LT.JLAST)  GO  TO  216 

212  MUPEX=MURcX 

GO  TO  (213.214.215) .MUREX 

213  NGEOM= 1 

CALL  GEOM 

DUDT ( J ) *PQ ( JMH ) /GRAMS ( JMH ) /2,  *AREA 
GO  TO  220 

214  DUDT ( J ) =0. 

C  GACCEL  TERM  CAN  CAUSE  MOTION  OF  LAST  INTERFACE  EVEN  WHEN  DUDT«0 

C  FIX  THIS  BY  TEMPORARILY  SETTING  GACCEL-O. 


HYDR 
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GSAV£=GACCEL 

GACCEL=0. 

GO  TO  220 

HYDR 

215 

CALL  BDY2 

GO  TO  220 

216 

NGEOM= 1 

CALL  GEOM 

DUDT ( J ) = ( PQ ( JMH ) -PQ ( JPH )  )  /  < GRAMS ( JMH ) +GRAMS ( JPH )  ) *2 • *AREA-GACCEL 
C  ADVANCE  U  FROM  TIME  N-l/2  TO  TIME  N+l/2 


220  U< J)=U< J)+DTMIN(3)*DUDT< J) 

C  ADVANCE  X  FROM  TIME  N  TO  TIME  N  +  l 

X< J)=X< J)+DTMIN(2)*U( J) 

NGEOM=2 

CALL  GEOM 
PWAS I JMH ) =P ( JMH ) 

EWAS< JMH)=E( JMH) 

DWASI JMH)=D( JMH) 

C  ADVANCE  D  FROM  TIME  N  TO  TIME  N+l 

D ( JMH ) =GRAMS ( JMH ) /VOL 

223  I F ( UCUT .LE.O* )  GO  TO  226 

C  SET  VELOCITIES  LESS  THAN  UCUT  EQUAL  TO  ZERO 

224  IF(ABS<U(J)).LT.UCUT)  U(J)=0. 

C  NDOALL= 1  FORCES  CALCULATION  OF  EVERY  20NE  IN  THE  PROBLEM. 

226  I F { NDOALL. EQ. 1 )  GO  TO  231 

227  IF( J.LT.JCALC-4)  GO  TO  231 

228  IF(U< J) .NE.O. )  KU=~1 

229  KU=KU+1 

CALCULATE  ARTIFICIAL  V I 5COS I TY ♦  Q. 

SAVE  PREVIOUS  VALUE  OF  0 

231  QWAS=0(JMH) 

IF(ABS( (D( JMH)-DWAS(JMH) )/DWAS(JMH> ) .GT.2.E-8)  GO  TO  2320 
D( JMH)=DWAS< JMH) 

GO  TO  238 

2320  IF(D< JMH) .LE.DWASt JMH) )  GO  TO  238 

232  VOLSUM=l ./D( JMH)+1./DWAS( JMH) 

X0IF=X( J)-X( J-l ) 

2321  IF(NC(7).GT.O)  GO  TO  2323 

C  USUAL  FORM  OF  VELOCITY  GRADIENT. 

2322  UGRAD=U< J-l ) -U ( J ) 

IFlUGRAD.LT. 0. )  GO  TO  238 
GO  TO  2324 

C  ALTERNATE  FORM  OF  VELOCITY  GRADIENT. 

2323  UGRAD=XD1F/DTMIN( 3)*( 1,-DWASI JMHJ/DI JMH) ) 

2324  IF(CINOd)  .LE.O. )  GO  TO  234 

233  Q( JMH) =2.*CINQ< I ) **2*UGRAD**2 /VOLSUM 
IF(AINQ( I ) .LE.O. )  GO  TO  239 

GO  TO  235 

234  Q ( JMH ) =0 • 
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235  QL I  N  =  A I  NO ( I ) «UGRAD/ VOL SUM«SOUND ( JMH )*2 . 

2 37  0( JMH)=0( JMHJ+QLIN 
GO  TO  239 

238  0 ( JMH 1=0* 

C  00  NOT  CALCULATE  MORE  ZONES  ?F  THE  LAST  4  ZONES  HAD  ZERO  U. 

C  KU  COUNTS  THE  CONSECUTIVE  ZONES  OF  ZERO  U. 

239  NOST  =NEQST ( I ) 

CALL  EOST 

PO ( JMH ) =P ( JMH ) +0  <  JMH ) 

T F ( JUTEST.GE. J)  GO  TO  240 

1 F ( NDOAL  L«  EO. 1 )  GO  TO  240 

IFIKU.LT.4)  GO  TO  240 

KU  =  0 

JCALC=J 

JCALC1= J-l 

JMAX=J 

GO  TO  251 

240  CONTINUE 

250  CONTINUE 

END  MAIN  HYDRODYNAMIC  LOOP 

NEXT  STATEMENT  IS  REACHED  ONLY  WHEN  ALL  ZONES  WERE  CALC. 

I F ( NDOALL . EQ. 1  )  GO  TO  251 
JCALC=JLAST 
JCALC1=JL AST-1 

RESET  GACCEL  TO  ITS  PROPER  VALUE 

251  GACCEL=GSAVE 

BEGIN  TIMESTEP  CALCULATION 

260  DO  280  1  =  1  *  1  MAX 
1=1 

JM I N= I NTER J ( I )+l 
JMAX  =  INTER J  ( l+l ) 

DO  275  J  =  JM I N « JMAX 
JMH=J-1 

IF( J.GT.JCALC)  GO  TO  201 
262  UD I F=ABS  (U( J) -U< J-l ) ) 

DENOM=SOUND( JMH)  +4 . *C I NQ ( I ) ** 2*UD I F  +2.*A INQ ( I ) #SOUND ( JMH ) 
STABIL  IS  A  STABILITY  CONSTANT  (ABOUT  .8). 

DTZ J ( JMH ) =STAB I L* ( X ( J )-X ( J-l ) ) /DENOM 
275  CONTINUE 

280  CONTINUE 

FIND  ZONE  WITH  SMALLEST  TIMESTEP. 

281  JDT=  1 
DO  283  J  =  2  » JCALC 1 
IF(DTZJ(J).LT.DTZJ( JDT) )  JDT*J 

283  CONTINUE 
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DTN£W=DTZJ< JOT) 

C  SAVE  PREVIOUS  VALUE  OF  T ! ME STEP 

28A  DTM1N(1)=0TM!N<21 

NEW  TtMESTEP  WILL  BE  USEO  TO  ADVANCE  X. 

DO  NOT  LET  0TMINI21  EXCEtU  DTRATE  TIMES  PREVIOUS  T 1MESTEP • 
OTMIN(2)=AMINl(DTNEW,OTRATt*DTMIN{ 1  )  ) 

AVERAGE  OF  OLD  ANu  NEW  7 IMLSTEPS  WILL  BE  USED  TO  ADVANCE  U. 
DTMIN(3)=(OTMIN(2)+DTMIN( I) )/2. 

END  TIMESTEP  CALCULATION 

299  RETURN 
END 


HYDR 
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tIBFTC  OUT  1  LIST  .REF 
SUBROUTINE  0UT1 


0UT1 


####*  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 
C 

DATA  NOUGHT .MCY.NTYMES/O .0 .0/ 

C 

I F ( IDENT.NE.l >  GO  TO  324 

301  WR I TE ( 6  » 302  ) 

302  FORMAT ( 40H0SUBRQUT 1 NE  0UT1  VERSION  3,  3-  1-65  ) 

ZNC6-NC ( 6 ) 

TFACTR=10.**! 1./ZNC6) 

N0UGHT=0 

RETURN 

C 

303  FORMAT ( 1H1 * 17HPROBLEM  NOL  WUNDY.I5.2H  ,12A6> 

304  FORMAT ( I5.1P6E13.5.2E8.1.E9.2.E12.5) 

305  FORMAT ! I6.1P3E14.5.I7) 


306 

FORMAT! 119H0 

J  X ! J+ 1 J 

U !  J  +  l )  X 

CENTER 

DENSITY 

ENERGY 

1  PG(J) 

0 ( J !  TOTPSI 

OR  ZONE  MASS 

TEMP 

DTZJ 

GAMMA J 

307 

FORMAT! 119H 

CM 

CM/SEC 

CM 

GM/CC 

ERG/GM 

1  DYNE/SOCM 

DYNE/SOCM  OVERPSI  GRAMS 

KELVIN 

SEC 

) 

308 

FORMAT! 119H0 

J  X ( J+ 1 ) 

U !  J  +  l )  X 

CENTER 

DENSITY 

ENERGY 

1  PO ! J ) 

Q  « J I  TOTPSI 

OR  DYNPSI  OR 

TEMP 

DTZJ 

GAMMA) 

309 

FORMAT! 119H 

CM 

CM/SEC 

CM 

GM/CC 

ERG/GM 

1  DYNE/SOCM 

DYNE/SOCM  OVERPSI  G/SO  CM 

KELVIN 

SEC 

) 

312  FORMAT! 78H0REGIOM  MATERIAL  NZONES  GAMMA  K-ENERGY 

1  I -ENERGY  TOT  ENE  ) 

313  FORMAT! I  5 . 2  I  8 . 4E 1 5 . 5 ) 

314  FORMAT! I6.E15.5) 

315  FORMA T ! 55H  CYCLE  TIME!SEC)  NEXT  TIMESTEP  TOTAL  ENERGY  JDT ) 


316  FORMAT! 50H1THE  FOLLOWING  DATA  HAVE  BEEN  WRITTEN  ON  TAPE  18-  J 

317  FORMAT! 53H1  *****  ENERGY  CHECK  *****) 

318  FORMAT! 1 4 , 1 P 1 E 1 2 . 5 , 1P5E 1 0 . 3 , 1 P 1E9 . 2 , 1 P3E 1 0 . 3 . 1 P 1 E9. 2 . OPF 5 . 2  > 

319  FORMAT ! 1 1 6H0  J  X  U  D  E 

1  0  PQ  TEMP  DTZJ  ) 

320  FORMAT! 14, 1P8E14.5) 

321  FORMAT! 16) 


GO  DIRECTLY  TO  PRINTING  IF  TROUBLE  HAS  OCCURRED. 
324  IF ( KUTOFF.NE.O )  GO  TO  346 


TAPE  18  IS  WRITTEN  HERE  IF- 

1)  NCYCLE=0.  OR 

2)  NTAPE  CYCLES  HAVE  PASSED  SINCE  LAST  WRITING.  OR 

3)  NWRITE  IS  POSITIVE. 

I  F ( NTAPE . LE .0  I  GO  TO  335 
326  NTYMES=NTYMES+1 

IFINCYCLE.EO.O )  GO  TO  329 
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C  NWRITE  POSITIVE  FORCES  WRITE  TAPE  18  UNLESS  NTAPE=0.  OUT1 

I F ( NWR I T£«GT • 0  )  00  TO  329 

C  NTYMES  COUNTS  THE  CYCLES  SINCE  THE  LAST  TAPE  WRITE, 

IF(NTYMES.LT.NTAPE)  GO  TO  33S 
329  NTYMES=0 

NRECRD=NRECR0+1 

JLAST1=JLAST-1 

C 

WRITE! 18)  T,(DTMIN(L) ,L=1 ,3) , NCYCLE ,NRECRU, JLAST , JCALC  ,  JCALC1  »*  . 


1  K.REZ3B ,  I  MAX  » NPROB  * 

WRITE! 18 )  (NZONES! I ) » INTER J! 1  +  1 ) . 1=1, 1  MAX ) ,X( JLAST  > ,U( JLAST ) *  * 

1!X(L) »U  ( L  ) , D  t  L ) »E  ( L )  »Q( L ) *PQ( L ) . TEMPI L) *DTZ J ( L ), L= 1 » JLAST 1 )  * 

NWR I TE=0 

IF(NC( 19) .LT.l  )  GO  TO  335 

WRITE16.316)  * 

WR I T£ ( 6  » 330  )  * 

330  FORMAT  !?3H0  T  NCYCLE  NRECRD  JLAST  JCALC  JCALC1  K.KE 

1Z3B  IMAX  NPROB  ) 

WR  I  T  E  (  6  *  332 )  T , NCYCLE, NRECRO, JLAST  * JCALC » JCALC 1 »<REZ3B  » l MAX  *  * 

1  NPROB  * 

332  FORMAT! 1H0,1P1E14. 5, 817) 

WRI TE (6*331 )  ( D  T  M I N ! L ) » L  =  1 »  3  )  * 

331  FORMAT! 1 OHODTM I N ( L ) = 1 P3E 16 . 7 ) 

WR I T  E ! 6 i 3  3  3 )  ( NZONES ( I ) ♦ I = 1 » IMAX )  * 

333  FORMAT! 11HONZONES! I )=20I5) 

I  MAX  1 = 1 MAX+ 1 

WRITE(6,334)  (INTER J(I)»I  =  1» I  MAXI )  * 

334  FORMAT ( 1 1 H  I NTER J < I  ) = 20  I  5 ) 

WRITE(6,319)  * 


WR I T  E ( 6 , 320 )  (L.X(L),U(L)*D(L)*E(L>,0(L),PU(L) » TEMP ( L ) ,DTZJ(L ) »  * 

1  L  =  1 ♦ JL AS  T 1  )  * 

WR I TE ( 6 ,320  )  JLAST , X < JLAST ), U ( JLAST )  * 

C 

335  I F ( NPR I  NT • GT • 0 )  GO  TO  341 
I F ( NCYCLE.GT  »  5 )  GO  TO  338 

336  MCY=MCY+1 
GO  TO  341 

C 

C  PRINTOUT  OCCURS  HERE  IF  LITHER- 

C  1)  NPR  CYCLES  wERt  CALCULATtO  SINCE  THE  LAST  PRINTOUT,  OR 

C  2)  SOME  SUBROUTINE  HAS  SET  NPRINT=1,  OR 

C  3)  NCYCLE  IS  LESS  THAN  5,  OR 

C  4)  TIME  T L I S T ! 7 )  HAS  ULLN  RLA^nED,  OR 

C  5)  TROUBLE  HAS  OC^JRREu  < KUTOFF = 1 ) , 

C 

338  IF(T.GE.TLIST{7) )  GO  TO  3400 

C  MCY  COUNTS  THE  CYCLES  SINCE  THE  LAST  PRINT 

339  MCY=MCY+1 

IF (MCY. LT. NPR )  RETURN 
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3400  2 F ( NC ( 6 ) )  340*340.3410 
3410  TLIST(7)*TLIST ( 7 ) *TFACTR 
MCY*0 
60  TO  346 

340  MCY=0 

341  NPR I N  T  =  0 

QUTDT  =  TL I  ST  < 1 ) 

1F(TLIST(7).GE.TL1ST(2)  )  OUTDT  =  TL I  ST ( 3 ) 

345  TL 1ST ( 7 ) =T+OUTOT 

CALC  TOTAL  INTERNAL  AND  KINETIC  ENERGY  IN  EACH  REGION. 

346  EINSUM=0» 

EKSUM=0. 

ESUM*0. 

DO  355  1 = 1 . IMAX 
El  NT ( I ) »0. 

EK I N ( I ) =0* 

JM 1 N* I NTER J ( I  ) 

JMAX= INTER J ( I +1 ) -1 
DO  350  J  =  JMI N  *  JMAX 
JPH= J 

USQ=  <  U ( J ) **2  +Ut J+l)**2)*.5 
EINT(  I  )*EINT (  I  )  +  E<  JPH)*GRAMSUPH) 

EKIN(I)=EKIN( I )  + • 5*USQ*GRAMS ( JPH } 

PDYN ( JPH )  =  • 5*D ( JPH ) *USQ*14. 5Q382E-6 
350  CONTINUE 

EINSUM=EINSUM+EINT(  I ) 

355  EKSUM=EKSUM+EK I N ( I ) 

ETOTAL=EINSUM+EKSUM 

ETENTH=.1*ESTART 

SUPPRESS  ENERGY  CHECK  IF  NC(4)=1 
IF ( NC ( 4 ) «E0« 1 )  GO  TO  360 

357  IFIABS  (ESTART-ETOTAL)-ETENTH)  360.360.358 

TERMINATE  PROBLEM  IF  TOTAL  ENERGY  CHANGES  TEN  PERCENT. 

358  WRI TE ( 6  »3 17 ) 

KUTOFF= 1 


360  WRI TE l 6 * 303  )  NPROB . < AL I S T{ L ) . L= 1 3 . 24 ) 

363  WR I TE ( 6 » 3 1 5  ) 

WRITE (6, 305)  NCYCLE.T.DTMINI2) .ETOTAL.JDT 
C 

JPRINT=JCALC-1 

C  NC(17)*1  CAUSES  PRINTING  OF  CUMULATIVE  MASSES  (GRAMS/SQ  CM 

C  BETWEEN  ORIGIN  AND  X)  IN  PLACE  OF  PDYN. 

IF(NC(17).LE.0)  GO  TO  3639 
IF(NCYCLE.EO.O)  JPR INT= JLAST-1 
PDYN( 1 )  =  D  ( 1 ) *  f  X ( 2 ) —  X ( 1 )  ) 

DO  3631  JPH=2  * JPR INT 

3631  PDYN( JPH )  = PDYN (JPH- 1 ) +D < JPH)* (X( JPH+1 )-X( JPH) ) 
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GO  TO  366 

C  MAS  = 1  CAUSES  PRINTING  OF  ZONE  HA5SES  IN  PLACE  OF  POYN. 

3639  IF(MAS.LE.O)  GO  TO  369 
JPRINT=JLAST-1 

364  DO  365  JPH=l,JPRiNT 

365  PDYN ( JPH ) =GRAMS ( JPH ) 

366  WR I TE ( 6  *  306  > 

WR  l  TE ( 6 « 307 ) 

MAS  =  0 

GO  TO  370 
C 

369  WR I TE ( 6  » 308 ) 

WR I TE ( 6  *309 ) 

370  IF(NCYCLE.EQ.O)  JPRINT=JLAST-1 
DO  361  J  = 1 ♦ JPR I  NT 

XAV( J  >  =  { X { J)+X(J+1) )*.5 
361  CONTINUE 

DO  371  L=1 s JPR INT 

371  PSI (L)=PQ!L) *14.50382E-6 
IFINC(II))  377  ,377,3710 

C  CALC  OVERPRESSURES  IF  NC(11)=1  (3710  THRU  375). 

3710  DO  375  L=1,JPRINT 

I F ( I ARDC.GE.3l )  GO  TO  372 

ALTCM=C0SPHI*(X(L)+X(L-1) >*.5  +H0BKM*1.£5 
CALL  ARDC(ALTCM,PRATIO,TRATIO,DRATIO) 

PR AM8  =  PRAT  1 0*14 .696 
GO  TO  373 

372  PRAMB=PQ( JLAST-1 I *14. 50382E-6 

373  PSI (L)=PSI (L)-PRAMB 

I F ( ABS ( PS  I (L) )*1.E  +  5.LT .PRAMB  )  PS1 (L)=0. 

375  CONTINUE 

377  WR I TE ( 6 , 3 18  )  NOUGHT , X ( 1 ) ,U ( 1) 

JSTART=1 

DO  385  11=1, IMAX 

JQU I T=JST ART+N ZONES ( I  I ) -1 

I F ( JQUIT.LT. JPRINT)  GO  TO  380 

I F ( JSTART.GE.JPRINT )  GO  TO  386 

379  JQU I T  =  JPR I  NT 

380  WRITE(6,318)  ( L ,X ( L+l ) ,U ( L+ 1 ) , XAV ( L ) ,D ( L ) , E ( L ) * PQ ( L ) »Q ( L ) » 

1  PS  I  (L) »PDYN ( L ) ,TEMP(L) «DTZJ( L ) , GAMMA J( L) *L=J START  » JQUI T ) 

WR I TE ( 6  *  38 1  ) 

381  FORMAT! 1H  ) 

JSTAR  T  =  JOU I T+ 1 

385  CONTINUE 

386  WR I TE ( 6 , 3 1 2  ) 

DO  390  1=1, IMAX 

ENTOT ( I ) =EK I N ( l  ) +E I  NT ( I ) 

390  CONTINUE 

391  WRITE  (6,313)  (I  ,NEQST  (  I  )  ,NZONES(  I  )  ,  GAMMA  I  (  I  )  *  EK.  I N  (  I )  »  E I  NT  (  I ) » 

1  F.NTOT(  I), 1  =  1, IMAX) 


OUT  1 


« 


* 

* 


* 


* 

» 

« 


* 


# 


* 
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C  LOCATE  ZONE  JOMAX  WHICH  CONTAINS  MAx  Q«  OUTl 

QMAX=0. 

JQMAX=0 

DC  393  J=  1 » JCALC 

1F{Q< J) .LE.QMAX)  GO  TO  393 

392  QMAX=Q ( J I 
JQMAX=J 

393  CONTINUE 

IF ( JOMAX *£0.0 i  GO  TO  399 

394  IF ( JOMAX »GE.J LAST -2 )  GO  TO  399 

C  INTERPOLATE  WITHIN  ZONE  JOMAX  FOR  t>E  T  TER  0  LOCATION, 

396  QFRAC=(0( JOMAX )-w ( JOMAX- 1 ) ) / { 0 ( JOMAX 1 -0 ( JQMAX+ 1 ) + 1 . E-l 5 > 

XQMAX=<  X( JOMAX )+X{ JOMAX+1 ) *OFRAC ) /{OFRAC+I , ) 

397  WR I TE ( 6 , 398  )  JOMAX, XQMAX  * 

398  F  ORMAT ( 2  3H0MAX I  MUM  G  IS  IN  ZONE, 14, 6H  AT  X=*1PE11.4) 

C  FIND  MAXIMUM  PRESSURE. 

JPMAX  =  JOMAX 
DO  405  J=l, JCALC 

I F ( PO ( J  J ,GT.PQ{ JPMAX  ))  JPMAX= J 

405  CONTINUE 

4050  XPMAX=XAV( JPMAX ) 

WR I TE ( 6  * 406 )  JPMAX, XPMAX  * 

406  FORMAT ( 2 3H  MAXIMUM  PSI  IS  IN  ZONE, 14, 6H  AT  X*,1P£11.4) 

C  STOP  PROBLEM  IF  TROUBLE  HAS  OCCURRED. 

399  IF (KUTOFF.EO.O)  RETURN 

IF(NC(12 J.GT.O)  CALL  DUMP ( A I NQ(1 ) ,ZONS I Z ♦ 1 , A I  NO (1) ,ZONS I Z *  2  ) 

STOP 

END 
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S1BMC  OUI2  US!  .REF 
SUbKUUI INE  OUT 2 


OUT  2 


L 

L 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


USUALLY 

KOUI  2 

KOUT2 

KOUT2 

KOUT2 

KOUT2A 

KOUT2B 


KOUT2B 
NC  ( 9 ) 

NCI  18) 


,  <0UI2=4» 
=  0 
=  +l 
=  +  2 
GR+2 

=  L 


=  0 


KOU 1 2A= 10  *  K0UT2B  =  0 ♦  NC I  9 ) =Q .NC l 18 ) 31 • 

THIS  PROGRAM  NOT  USED. 

OMIT  tXTRAPOLAT I  ON  OF  PRESSURE. 

INCLUOt  tXTRAPOLATLO  PO  IN  OUTPUT. 

BASE  THE  EXTRAPOLATION  ON  K.0UT2  POINTS  PER  GROUP. 
PREPARE  A  LINE  OF  OUTPoT  tVERY  XOUT2A-TH  CYCLE. 

USE  S ( 1  THRU  L)  AS  OVERPRESSURES  AND  FIND  THE  DIST¬ 
ANCES  AT  WHICH  THEY  OCCUR.  WHEN  A  PRESSURE  SIN)  IS 
REACHED.  REPLACE  SIN)  BY  THE  DISTANCE  AND  INCREASE 
NCI26)  BY  1.  ITHESE  DISTANCES  ARE  BEING  FOUND  FOR 
USE  BY  OUT3). 

OMIT  THE  SIL)  SECTION. 

=  0  PUT  OUT 2  OUTPUT  ON  TAPE  6, 

=2  PUT  OUT 2  OUTPUT  ON  TAPE  26. 

DO  NOT  CARRY  MAX  G  SEARCH  INWARD  INTO  REGION  NC I  1 8 ) 


DIMENSION  MSY (50)*TAIM{50) .MJ I  50  S »  XPMAX I  50 ) »OVPQ I  300 ) . JBACXl 50 ) 
1  .CVPQEXI 50) .AMBPSI ( 50) .DLPDRI50) . XQMAX I  50 ) .OVPMAX ( 50 ) 


OUT  2 
OUT  2 


*****  INSERT  STANDARD  COMMON  CARDS  dEFQRE  COMPILING  ***** 

C 

DATA  NUFSED/0/ 

C 

IF  I IDENT.NE.l )  GO  TO  7100 

701  WRITEI6.702) 

702  FORMAT  1 50H0SUBROUT INE  OUT2  VERSION  3.  3-  1-65  TAPE  6.26  ) 

7100  IFINUFSED.E0.123A)  GO  TO  710 

NUFSED= 1 2  3A 


XMAXQ=0. 

QHUGE= 1 . 

I COUNT=0 
M00  =  0 
NP  T  S  =  4 
NC I  26 ) =0 

I F I KOUT 2-2  )  710.704 .705 

704  N I L E  =  4 

GO  TO  706 

705  NI LE  =  KOUT 2 

706  NILE2=2*NILE 
SN I L£=N I L£ 

RETURN 

710  I COUN T  = I  COUNT* 1 

I F I NOUI T.LT.O )  GO  TO  716 

715  IF( ICOUNT.LT. KOUT2A)  RETURN 

716  I COUN  T  =  0 
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C  LOCATE  ZONE  JQMAX  WITh  HIGHEST  0.  0UT2 

C 

QMAX=0 • 

JOMAX=] 

NC 18=NC ( 18  ) 

JBEGIN=INTERJ(NC18+1> 

JCALC1= JCALC-1 
DO  721  J=JBEGIN,JCALC1 
I F ( Q ( J ) *L£ .QMAX )  GO  TO  721 

720  CiMAX  =  Q(J) 

JQMAX= J 

721  CONTINUE 

INTERPOLATE  WITHIN  ZONE  JQMAX  FOR  bETTER  Q  LOCATION. 

USc  MIDDLE  OF  ZONE  IF  JQMAX  IS  AT  EDGE  OF  PROBLEM. 

IF (JQMAX. LE. 1 )  GO  TO  723 

722  IF( JQMAX. LT.JLAST-l )  GO  TO  726 

723  XMAXO= ( X ( JQMAX+1 ) +X ( JQMAX ) )*.5 
GO  TO  726 

726  QFRAC-IQ; JQMAX )-Q{ JQMAX- 1 ) )/( Q ( JQMAX ) -Q ( JQMAX+1 J+l.E-15) 

XQWAS=XMAXQ 

XMAXQ=( X( JQMAX )+X ( JQMAX+1 ) *QFRAC ) / ( QFRAC+1 • i 
1 F ( XMAXQ. GE . XQWAj )  GO  TO  726 
XMAXQ=XQWA5 
GO  732 
.26  MOO=MOO+ 1 

XCMAX ( MOO ) =  XMAXQ 
MSY(MOO)=NCYCLE 
TAIM(MOO) =T 
MJ ( MOO )= JQMAX 

CALCULATE  AMblENT  PRESSURE  AT  LOCATION  OF  MAX  Q. 

I F ( 1 ARDC-3 1 )  729,728,728 

728  AMBPSI (MOO)=PQ( JCALC-1 ) *  14 . 50 J 82 1-6 
GO  TO  7301 

729  ALTCM  =  COSrHI  +XQMAX  (  MOO  )  +HOPK.M*  1  •  L  6 
CALL  APDCI ALTCM,PRATIO.TRATIO,DRATIO) 

AMBPSI { MOO ) =I4.69b*PRAT  10 

SEARCH  FOR  LARGEST  OVERPRESSURE. 

C  CALCULATE  OVERPRESSURES.  (THRU  7305). 

7301  DO  7305  L  =  J3EG  IN. JCALC1 

I F ( I ARDC-3 1  )  7303 ,7302, 7302 

7302  AMBPR=PO( JCALC-1 1+14.50 382E-6 
GO  TO  7304 

730  3  ALTCM-COSPHi  »  (  X  (  L  +  l  )+X(L)  )  *  •  5  +HOBK.M*  1 ,  E  5 
CALL  ARDC (ALTCM.PRAT lO.TRAT iO.DRAI 10) 

AMBPR= 14.696+PRAT 10 

7304  0VPQ(L)=PQ(L)»14.50382E-6  -AMbPR 

7305  CONTINUE 
C 
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730  PGBIG=OVPQ< JBEGIN)  0UT2 

JPQMAX= J8EG I N 
DO  735  J=  1  * JCAIC1 
IF(0VPQ< J) .LE.PQBIG)  GO  TO  735 
POBIG=OVPO( J) 

JPQMAX=„ 

735  CONTINUE 

JSAVE=JQM AX- JPQMAX 

MAXIMUM  OVERPRESSURE  WAS  FOUND  IN  ZONE  JPQMAX. 

JSAVE  IS  NUMBER  OF  ZONES  BACK  FROM  MAX  Q  TO  MAX  OVPQ. 

IF (JPQMAX. LE.JBEGIN)  GO  TO  732 
IF( JSAVE. LT. 24)  GO  TO  739 
732  XPMAX(MOO) =-0. 

OVPMAX(MOU)  =  -0. 

OVPOE  X ( MOO ! =~0  « 

JBACK! MOO) =0 
GO  TO  /005 

739  X PMAX ( MOO )  =  ( X ( JPQMAX ) +  X ( JPQMAX+ 1 )  )*.5 

740  OVPMAX(MOO)  =  POBIG 
JBACK ( MOO ) = JSAVE 

TEST  IF  EXTRAPOLATION  UESIREO 

IF(KOUT2.LT,2)  GO  TO  753 

CALCULATE  EXTRAPOLATED  OVtRPRESSURE  oY  EXTRAPOLATING  P  VS  X  DATA 
TO  THE  POSITION  OF  MAXIMUM  U. 

741  I  F  (  JPQMAX- JbEG  IN.UT  . N  I  Lt<;  >  UO  TO  743 

742  OVPOEX ( MOO )  =  -0. 

GO  10  753 

743  SUMA=0. 

SUMb-0 » 

DO  745  J= 1 *N  I  L  E 
J J  = JPQMAX+ 1 - J 
$UMA=SUMA+OVPQ( JJ ) 

SuMt>  =  6UMB  +  X  (  JJ  ) 

745  CONT INUE 
SUMC=0. 

SUMU=0 . 

DO  748  J= 1 ♦ N I L  E 
JJ=JPOMAX-NILE+l-J 
SUMC  =  SUMC  +OVPQ ( J J ) 

SUMU  =  SUMD+X ( JJ  ) 

748  CON  T INUt 

PA VI =SUMA/SN I Lh 
XAVl=buMB/sNlLE 
PAV2=bUMC/SNILE 
XAV2=SUMD/SNILE 
SLOPE = ( PA VI -PAV2 ) / <  XAV 1-XAV2 ) 
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OUT  2 


TERCEP=PAV2-SLOPE«XAV2 
OVPOEX  (  MGO  )=  SLOPE* XUi'.AX!  MOO  )+T£R CLP 

SPECIAL  OPTION  FJR  TRIGGcRING  OUT 3  AT  uESIRtD  PRESSURE  LEVELS-. 

K.OUT 2b=NuMbEk  OF  PRESSURc  LEVELS.  RtAD  IN  ON  S  CARDS. 

753  I F ! KOUT28 )  7005,7005,7000 

7000  NLOC  =  NC !  26 )  +1 
IF (OVPOEX (MOO) )  7005*7005,7001 

7001  IF (OVPOEX (MOO) .GT.S(NLOC) )  GO  TO  7003 

S ( NLOC+1 2 ) = AMO PS  I  ( MOO )  C 

7002  S ( NLOC ) =XuKAX ( MGu ) 

NC ( 26 ) =  NLOC 

I  F  (  K.0UT 2B.LE.NLOC  )  KOuT2b  =  0 

PRINT  WHEN  50  LINES  OF  DATA  HAVE  btEN  STORED  OR  IF  NQUIT  IS 
NEGATIVE  (INDICATING  COMPLETION  OF  PROGRAM) 


7005  IF (NOUI T.LT.O)  GO  TO  761 
756  I F ( MOO j  L  T «  50 )  RETURN 

761  I P { NC  <  9 ) .GT .0 )  GO  TO  7610 

WRITE!  6,762)  * 

WR I TE (  6,7620)  * 

WRITE!  6,763) 

GO  TO  7 6<* 

7610  WRITE! 26,762)  * 

WRI TE(26,7620)  * 

WR  I  TE  (  <!6 , 763  )  * 

762  FORMAT! 17H1SHOCK  FRONT  DATA  ) 

7620  FORMAT! 1HO,31X,A2HcXTkAPOLaTcD  VALuES  DIRlCTLY  READ  VALUES  ) 

763  FORMAT! 119H  CYCLE  J  JbACX  TIME  XIMAX  0)  MAX  OVP  X! 

1MAX  OVP)  MAX  OVP  SLCPc  AMbPSI  > 

MOO=NUM3ER  OF  LINES  OF  DATA  STORED. 

764  DO  790  J=1 ,MOO 
IFtJ.LT. 5)  GO  TO  781 

IF ( J.GT.MOO-5 )  GO  To  781 


CALCULATE  LOGARITHMIC  SLUPc  OF  gVcKPkLSSURE  VS  DISTANCE  CURVE 
USING  4  POINTS  ON  EACH  SlUt-  OF  ulSIRtu  POINT. 

766  PSUM 1=0. 

PSUM2  =  0 • 

XSUM 1=0. 

XSUM2=0. 

DO  780  J J= 1  * 4 
JM= J~JJ 
JP=J+JJ 

IF  (OVPMAX(JM)  .LE.  0.)  GO  TO  781 
IF  (OVPMAX(JP)  .LE.  0.)  GO  TO  781 
779  PSUM1  =  PSUM1  +  ALOGIO ( OVPMAX ( JM ) ) 

PSUM2  =  PSUM2  +  ALOG10(OVPMAX ( JP l ) 
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C 


XSUMl=XSUMl+ALOGlO ( XQMAX ( JM ) ) 
XSUM2  =  XSUM2  +  AI_OG10  <  XOMAX ( JP ) ) 

780  CONTINUE 
DL0GP=P5UM1-PSUM2 
0L0GX=XSUM1-XSUM2 
OLPDX=ULOGP/DLOGX 
GO  TO  782 

781  0LP0X=0 • 

782  DLPDR ( J ) =DLPDX 


WRI  TEl^TS^MSY  <  J)  *  M  J  (  J  )  » JtiACX  (  J  }  *  TA 1M( J )  ♦XUMAX(J)  »UVPQEX(J)  * 
1  XPMAX ( J ) ♦ OVPMAX ( J ) *DLPDR ( J ) ♦AMoPSl ( J) 

788  WRITE! 26«789 )  MSY ( J ) »MJ ( J )  ♦ JbACK ( J) »  TAIMI J  ) * XOMAX ( J )  iGVPQEX ( J ) * 

1  XPMAXI J) *OVPMAX( J) tOLPUKl J) *AMoPSI ( J) 

789  FORMAT  (  1H  ♦I6tIbvI<*,lP2E12.3*lPEl0«3*lPElA«6»lPE10»3t0PF7.3» 

1  1PE1A.5) 

790  CONTINUE 
M00  =  0 

799  RETURN 
END 


OUT  2 


* 

* 

* 

* 
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SIBFTC  OUT 3  LIST, REF 
SUBROUTINE  OUT3 


OUT3 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  SUBROUTINE  PREPARES  SUMMARY  SHhETS  OF  P,  U,  DENSITY,  AND 
DYNAMIC  PRESSURE  VERSUS  TIME  AT  UP  TO  8  DISTANCES,  S(l-8), 

IF  KOUT3=0»  THIS  SUBROUTINE  IS  NOT  CALLED. 

USUALLY,  NC ( 10  >  =0 »  NC{15)=0»  KOUT3  =  0  OR  4,  KOUT3A*l0  . 

NCI  10)  =0  PUT  OUT 3  OUTPUT  ON  TAPE  6. 

=2  PUT  OUT 3  OUTPUT  ON  TAPE  26. 

NCI  15)  =WHEN  SHOCK  ARRIVES*  STORE  DATA  EVERY  CYCLE  INSTEAD  OF 

EVERY  K0UT3A  CYCLES  FOR  NCI  15)  CYCLES.  USE  150  IF  0  . 

K0UT3  =  NUMBER  OF  POINTS  AT  WHICH  P  VS  T  IS  DESIRED  (FROM  0  TO  8) 

K0UT3A=  STORE  A  LINE  OF  DATA  EVERY  K0UT3A-TH  CYCLE. 

IF  K0UT2B  IS  NEGATIVE, 

THIS  SUBROUTINE  COMPUTES  THE  MOTION  (RADIAL  DISTANCE  FROM  BURST 
CENTER  VS  TIME)  OF  THE  POINTS  AT  WHICH  PT  DATA  ARE  BEING  TAKEN. 

AT  THE  BEGINNING  OF  THE  PROBLEM,  THIS  SUBROUTINE  EXPECTS  TO  READ 
DATA  CARDS  FOR  RA ,  Rb,  AND  SPEED.  THESE  CARDS  ARE  TO  BE  PLACED 
AFTER  THE  REGULAR  INPUT. 

RY  =  DISTANCE  (CM)  TO  NEAREST  POINT  ON  (STRAIGHT-LINE)  TRAJECTORY. 

RX  =  DISTANCE  I  CM )  ALONG  TRAJECTORY  MEASURED  FROM  NEAREST-TO-BURST 

POINT  IN  DIRECTION  OF  MOTION.  (I.E.,  ONCOMING  OBJECT  HAS 
NtGATIVE  RX,  OBJECT  GOING  AWAY  HAS  POSITIVE  RX.) 

SPEED  =  SPEED  I  CM/ StC )  OF  OBJECT  ALONG  TRAJECTORY. 

DIMENSION  RY( 8 ) ,RX( 8) , SPEED (8 ) ,RX IN  I T I  8 ) ,DIREC ( 160 ) »RL 1 160 )  0UT3 

1,  POPS  I  I  160) ,DENSTY( 160) *  TA i M ( 160  > , R I  8 ) ,MM ( 8 ) , NC YKLE 1 1 60 )  0UT3 

2. AM0PSI (8) ,0VPSI ( 160) ,UPT( 160) ,MST0RE(8) , MCOUNT ( 8 ) »NUFL ( 8 )  OUT 3 


****  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 


DATA  NONCE/O/ 

IF!  IDENT.NE.1 )  GO  TO  804 

801  WRITE(6,8P2) 

802  FORMAT ( 55HOSUBROUT I NE  OUT3  VERSION  3,  3-  1-65  TAPE  6,26  ) 

RETURN 


THIS  SECTION  (804  THRU  8)  IS  UStD  ONLY  ONCE. 

804  I F ( NONCE-1 2 34 )  805,807,805 

805  DO  806  LOC=l ,KOUT3 
MM  I LOC ) =0 

NUPH=  150 

I F ( NC ( 1 5 ) » GT .0 )  NUPH  =  NC (15) 

NUFL ( LOC ) =NUPH 
MSTORE(LOC) =KOUT3A 
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MCOUNT ( LOC ) *0  0UT3 

806  CONTINUE 

NL I N£S= 160/KOUT3 
I F ( NL INE5.GT • 50 )  NLINES*50 
NONCE- 1234 
KPT=KOUT3 
NPT=KOUT3 
I F ( KOUT2B )  1.9.899 
9  I F ( IARDC.GT.30)  GO  TO  6 
DO  10  L  =  1 .KOUT3 
R(L)«S(U 

ALTCM=COSPHI *R ( L ) +HOBKM* 1 . E5 

CALL  ARDCIALTCM.PRATIO.TRATIO.DRATIO) 


AMBPSI (L)=14.696*PRAT10 
10  CONTINUE 
GO  TO  807 

1  READ! 5.2)  (RY I L >  ,L  =  1 ,NPT )  * 

READ (5.25  (RX INI T ( L) ,L=1 »NPT )  * 

READ( 5.2)  ( SPEED ( L ) . L= 1 *NPT )  * 

2  FORMAT ( 1P6E10.5 ) 

I F ( NC ( 10) .NE.2 )  GO  TO  7 

WRITEI26.3)  (RY(L) ,L=1,NPT)  * 

WR  I  TE ( 26 ,4 )  (RXINIT(L) ,L  =  1.NPT)  * 

WR I TE ( 26  *5 )  ( SPEED ( L ) »L= 1 »NPT )  * 

GO  TO  6 

7  WR I TE ( 6 , 3 )  ;RY(L).L=1.NPT)  * 

WRITEI6.4)  (RXINIT(L) .L-l.NPT)  * 

WR  I TE ( 6  ♦  5 )  ( SPEED (L ) .L=l *NPT  J  # 

3  FORMAT ( 7H  RY  =1P8E14„5) 

4  FORMAT ( 7H  RX  =1P6E14.5) 

5  FORMAT ( 7H  SPEED=1P8E14. 5 ) 

6  DO  8  L=  1  *8 


AMBF>S1  (L)  =PQ(  JCALC  -1 )  *14.503  82E-6 
R ( L) =5( L ) 

8  CONTINUE 
GO  TO  807 

DATA  IS  STORED  FOR  LOCATION  LOC  EVERY  MSTORE(LOC)  CYCLES. 

MCOUNT ( LOC )=NUMBER  OF  CYCLES  SINCt  LAST  PRINT  FOR  LOCATION  LOC. 

807  I  F  (  K.0UT2B  )  812.812.808 

C  BYPASS  OUT 3  UNTIL  0UT2  HAS  SET  NCI26)  NON-ZERO. 

808  K.P T  =NC  (  26  ) 

IF(KPT.LE.O)  RETURN 

810  IFINUFL(LOC) .EO.NUPH)  MSTORE ( LOC ) = 1 
812  DO  890  LOC= 1 .KPT 

MCOUN  T ( LOC ) =MCOUNT ( LOC )  +  1 
IFINQUIT.LT.O)  GO  TO  815 
IFINCYCLE.LT. 5 )  GO  TO  890 

814  I F ( MSTORE ( LOC ) ,GT .MCOUNT { LOC ) )  GO  TO  890 

815  MCOUNT (LOC) =0 

C  SEEK  DATA  TO  PUT  IN  LINE  KK. 

KK=MM ( LOC ) + 1 
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KKSLOONL I  NES*  (  LOC-1 )  +KK 
I F ( KQUT28 )  816  ,833,832 
816  RX ( LOC ) =RX I N I T ( LOC ) +SPEED ( LOC ) *T 
RRR*SQRT  ( RY ( LOC ) **2+RX ( LOC )  **2 ) 

R(L0C)=AMIN1<RRR,X( JLAST-1)) 

GO  TO  833 

832  R ( LOC ) sS { LOC ) 

AMBPSI ( LOC j  *S  ( LOC+1 2 ) 

FIND  ZONE  JJJ  IN  WHICH  OESIRED  R  LIES. 

833  JJJ=JLAST-1 

DO  835  L  =  2  * JLAST 
IF(X(L)»LT»R(LOC) )  GO  TO  835 

834  JJJ=L-1 
GO  TO  836 

835  CONTINUE 

836  PHRAC* ( R ( LOC )-X ( J J J ) ) / ( X  C JJJ+1 )-X( JJJ) ) 

UPT ( KKSLOC ) =  U l JJJ )  +  PHRAC*(U( JJJ+1 )-U( JJJ ) ) 

NABOR  IS  THE  NEARER  OF  THE  TWO  ZONES  ADJACENT  TO  JJJ. 
IF{PHRAC.LT.0.5)  GO  TO  838 

837  NABOR *=  JJJ+ 1 
GO  TO  839 

USE  ZONE  2  AS  NABOR  IF  JJJ=1  . 

838  IF(JJJ.LE.l)  GO  TO  837 
8381  NABOR= JJJ-1 

839  CTR=(X( JJJ+1)+X(JJJ)  )/2. 

CTRADJ= ( X ( NABOR )+X ( NABOR+1 ) )/2. 

FRAC  IS  INTERPOL  FACTOR  FOR  UUANTITItS  DEFINED  AT  ZONE  CENTERS 
FRAC* ( R ( LOC ) -CTR ) / ( CTRADJ-CTR ) 

POPS  I ( KKSLOC )  =  { PQ ( JJJ ) +FRAC* ( PQ ( NABOR ) -PQ ( JJJ ) ) ) *  1 4 . 50382E-6 
RL ( KKSLOC ) *R ( LOC ) 

I F ( KOUT2B )  8390.8400.8400 

CALC  NET  PARTICLE  VELOCITY  INCLUDING  THAT  DUE  TO  STA  MOTION. 

8390  UTRA J=SPEED { LOC ) -UPT ( KKSLOC ) *RX ( LOC ) / R ( LOC ) 

UNORM=UPT ( KKSLOC ) *R Y ( LOC ) /R (LOC ! 

UPT ( KKSLOC ) =SQRT  (UTRAJ**2+UNORM**2 ) 

CALC  DIRECTION  COSINE  OF  FLOW  WITH  RESPECT  TO  TRAJECTORY. 

DIREC  IS  COSINE  OF  ANGLE  BETWEEN  STATION  DIRECTION  OF  MOTION  AND 
DIRECTION  OF  PARTICLE  MOTION. 

D I REC= 1  IS  HEAD-ON,  DIREC=0  IS  SIDE-ON,  DIREC=-1  IS  TAIL-ON. 
DIREC(KKSLOC)»UTRAJ/UPT (KKSLOC ) 

NOTE  THAT  THE  TRAJECTORY  IN  THIS  FORMULATION  REMAINS  A  STRAIGHT 
LINE  THROUGHOUT  THE  PROBLEM.  IF  DESIRED,  THt  TRAJECTORY  COULD  BE 
MADE  TO  VARY  ACCORDING  TO  CALCULATED  FORCES  ON  THE  STATION. 

8400  DENJJJ=D( JJJ) 


t 
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DENNAB=D ( NABOR )  OUTS 

DENSTY ( KKSLOC ) «DENJJJ+FRAC* ( DENNAb-DENJJJ ) 

TA I M  (  KKSLOC  )  =  T 

OVPSI (KKSLOC) =PQPSI ( KKSLOC ) -AMBPS I ( LOC ) 

I F ( ABS  (OVPSI { KKSLOC) )-.0001*AM8P5I (LOC) )  840.840*841 

840  OVPSI (KKSLOC) =0. 

841  NCYKLE ( KKSLOC ) =NCYCLE 

PRINT  WHEN  NLINES  LINES  OF  DATA  ARE  STORED  FOR  A  GIVEN  R. 

NOUIT  NEGATIVE  MEANS  PROBLEM  IS  DONE  AND  EVERYTHING  IS  PRINTED* 

MM ( LOC )  IS  THE  NUMBER  OF  LINES  OF  DATA  STORED  SO  FAR  IN  LIST  LOC. 

MM ( LOC ) =MM ( LOC ) + 1 
C  NUFL(LOC)  IS  NUPH  UNTIL  THE  SHOCK  ARRIVES. 

NUF=NUFL ( LOC ) 

IF(NUF)  846.8450*842 

842  IF (NUF.LT.NUPH)  GO  TO  845 

843  IF (OVPSI (KKSLOC) )  889.889.844 

844  MSTORE ( LOC ) =1 

845  NUFL(LOC)=NUFL(LOC)~l 

I F ( KOUT2B )  846.8451,8451 
8451  RY ( LOC ) =5 ( LOC ) 

GO  TO  846 

8450  MSTORE ( LOC ) =KOUT3A 
NUFL( LOC ) =  -l 
C 

846  IF(NQUIT)  848,847,847 

847  IF(KK.LT, NLINES)  GO  TO  890 

C  PRINTING  IS  FORCED  WHEN  848  IS  REACHED. 

848  I F ( NC ( l 0 ) • NE. 2  )  GO  TO  850 


8500  WRI TE( 26.851 )  LOC  * 

WR  I  TE  (  26 , 852  )  AMBPSKLOC)  * 

WRITE(26,8521 )  SPEED(LOC)  * 

WR  I  TE ( 26 ,8522 )  RY(LOC)  * 

WR  I  TE ( 26 , 853 )  * 

GO  TO  8530 

850  WRITEt  6,851)  LOC  * 

WRITE!  6,852)  AMBPSKLOC)  « 

WRITEt  6,8521)  SPEED(LOC)  * 

WRITE!  6,8522)  RY(LOC)  * 

WRITE!  6,853)  * 

851  FORMAT ( 10H1  STATION, 12) 

852  FORMAT ( 2  5H  AMBIENT  PRESSURE ( PS  I )  =  , 1PE12 . 5 ) 

8521  FORMAT ( 25H  SPEED ( CM/SEC )  =»1PE12.5) 

8522  FORMAT ( 2 5H  PERI-BURST  DIST(CM)  =,1PE12.5) 

853  FORMAT ( 106H0  CYCLE  TIME  PSI  OVPSI  DYN  PSI 

1  DENSTY  PART  VEL  DIREC  COS  RADIAL  R  ) 

8530  DO  860  LM= 1 , KK 

LMSLOC=NLINES» (LOC-1 ) +LM 

DYNPSI=.5*DENSTY( LMSLOC )*UPT(LMSLOC ) **2*14.50382E-6 
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I F { NC ( 10 ) »NE»2 )  GO  TO  8540  0UT3 

8541  WRITE!  26.854)  NCYK.LE <  LMSLOC  )  . TA I M |  LMSLOC )  *PQPS I «  LMSLOC )  »  OVPSKLM* 
1SL0C) .DYNP5I *DENST Y ( LMSLOC ) *UPT { LMSLOC  }  .01 RbC { LMSLOC ) .RL(LMSLOC)  * 

GO  TO  360 

8540  WRITE*  6*854)  NCYKLE  {  LMSLOC  )»  TA1  M(  LMSLOC  )  »PQPS  1  (  LMSLOC )  »  OVPSKLM* 
1SL0C) .DYNPS I .DENSTYl LMSLOC ) .UPT( LMSLOC) .01 REC( LMSLOC ) »RL {LMSLOC)  * 

854  FORMAT { 1H  *  1 6 . 1P8E12 . 3 ) 

860  CONTINUE 

889  MM ( LOC ) 20 

C  END  PRINTING. 

890  CONTINUE 

C  END  STATION  (LOC)  LOOP. 

899  RETURN 
END 
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SIBFTC  OUT 4  LIST  »REF 

SUBROUTINE  0UT4 

C  PUNCH  CAROS  WHEN  CYCLE  K0UT4  IS  REACHED.  0UT4 

C 

*****  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 

C 

IFUDENT.NE.l)  GO  TO  90 
WR I TE ( 6  »89 ) 

89  FORMAT (50H0SUBROUTINE  0UT4  VERSION  3,  3-  1-65  PUNCH  CARDS) 


RETURN 

C 

90  IFINCYCLE.NE.K0UT4)  RETURN 
KOUT4=0 
INC0DS=1 

PUNCH  1,  ( AL I  ST ( L ) *  L= 1 » 12 )  * 

PUNCH  1,  (ALIST(L) *L=13»24)  * 

1  FORMAT ( 1 2A6 ) 

PUNCH  2»  NPROB  »K  ,MURIN  ,MUREX  ,IMAX  ,IARDC  »  * 

1  I  NT APE » I NCODS  «NQU 1 T  ,NPR  »NT  APE  »JCALC  * 

PUNCH  2.  K0UT2  ,KOUT2A»KQUT2B,KOUT3  .K0UT3A.KQUT4  ,  * 

1  KREZ1  »KREZ2  ,KREZ3  .KREZ3A ,KREZ3B  ,KREZ4  * 

PUNCH  2.  (NC(L) »L=1 » 12 )  * 

PUNCH  2*  ( NC ( L ) »LS 1 3 » 24 )  * 

2  FORMAT { 1215) 


PUNCH  3.  (I  .NEOSTII)  .GAMMAI  II )  .NZQNESU)  .OUTBDYU)  .EINITI  I)  ,  * 

1  UINITI  I )  ,DINIT(  I  )  .CINOI  1  )  .AINQ(I)  .ZONINGU)  ,I«1,IMAX)* 


3  FORMAT ( 12, I3,0PF6.3,I4,1P4E10.3»0P3F5.2) 

PUNCH  4,  COSPHI .HOBKM.WKT  , BLANK 1 .BLANK2 » ZONS I Z  * 
PUNCH  4,  T,TREZ0,DTMIN(2) , DTRATE , ST AB I L »UCUT  * 
PUNCH  4,  ( TL I  ST ( I ), 1*1,6)  * 
PUNCH  4,  (S(L) »L= 1 » 12 )  * 


4  FORMAT ( 1P6E10.3) 

PUNCH  5 ,NC YCLE 

5  FORMAT ( 1 5H  SLIST1  NCYCLE= , 1 5 , 10H , I NC0DS  =  2S ) 
ZERO=0. 

JLAST1=JL AST-1 


PUNCH  91.  (L»X(L),U(L),D(L) » E  (  L  ) ,Q(L),L*1«JLAST1)  * 

PUNCH  91,  JLAST , X ( JLAST ) ,U(JLAST> .ZERO, ZERO, ZERO  * 

N300=300 

PUNCH  91,  N300, ZERO, ZERO. ZERO, ZERO, ZERO  * 

91  FORMAT! 1 4 , 1P4E 14. 7  * 1P1 E12. 5 ) 

NULL=0 

PUNCH  2.  NULL  * 

RETURN 
END 
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SIBFTC  REZ1  LIST »REF 
SUBROUTINE  REZ1 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

***** 

c 


REMOVE  SQUASHED  ZONE 

THIS  SUBROUTINE  CUTS  RUNNING  TIME  BY  COMBINING  THE  ZONE  WITH  THE 
SMALLEST  TIME  STEP  WITH  THE  ADJACENT  INNER  ZONE.  UNLESS  THIS 
WOULD  RESULT  IN  THE  REMOVAL  OF  AN  INTERFACE.  IN  WHICH  CASE  THE 
ZONE  IS  COMBINED  WITH  THE  ADJACENT  OUTER  ZONE. 

NOT  OCCUR  IN  REGIONS  OF  LESS  THAN  6  ZONES  OR  IN 
REGION  KREZ1.  KREZ1  SHOULD  BE  SET  SMALL  ENOUGH 
FINE  ZONES  OF  REZ3  ARE  INVOLVED.  (REMOVING  ONE 
ZONES  WOULD  FOUL  UP  REZ3 ) . 


RE21 


REZONING  WILL 
REGIONS  BEYOND 
SO  NONE  OF  THE 
OF  THE  FINE 


INSERT  STANDARD  COMMON  CARDS  bEFORE  COMPILING 


***** 


C 

C 


401 

402 


403 


404 


406 


I F ( IDENT.NE.15  GO  TO  403 
WR I TE ( 6.402 ) 

FORMAT ( 55H0SUBROUT I NE  REZ1  VERSION  3,  3-  1-65  TAPES  6.  26 

TBIG=0. 

RETURN 


TBIG  IS  LARGEST  TIMESTEP  USED  THUS  FAR 
1 F ( KREZ1 «EQ. 1 )  GO  TO  404 
IF(NCYCLE.EQ.KREZl)  GO  TO  412 
IF(DTMIN(3) .LE.T8IG)  GO  TO  406 
TB I G=DTMIN ( 3  I 
RETURN 

IF(DTMIN(3).GT.0.7*TBIG)  RETURN 


IN  PROBLEM, 


ZONE  JDT  HAS  SMALLEST  TIMESTEP 
DETERMINE  REGION  THAT  CONTAINS  ZONE  JDT. 

412  NJ= JDT 

DO  420  1=1.1 MAX 
IF»NJ-INTERJ( 1+1 ) )  418.419.420 

418  NJREG= I 
GO  TO  425 

419  NJREG=I +1 
GO  TO  425 

420  CONTINUE 

REGION! NJREG)  CONTAINS  ZONE(NJ),IF  REG  ION ( N JREG ) CONTAI NS 
LESS  THAN  SIX  ZONES. REZONING  IS  NOT  CARRIED  OUT 

425  IF<NZ0NES(NJREG).LT.6)  RETURN 

IF  NJREG  EXCEEDS  KREZ1.  REZONING  IS  NOT  CARRIED  OUT. 

IF(NJREG.GT.KREZ1 )  RETURN 
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TEST  TO  DETERMINE  WHETHER  ZONE(NJ)  IS  BOUNDED  A80VE  BY  A  REGION  REZ1 
INTERFACE*  S F  IT  IS  SET  NJ*NJ  +  1 

426  DO  427  I=1,IMAX 
IF(NJ.EQ.INT£RJ< I  1)  GO  TO  42S 

427  CONTINUE 
GO  TO  429 
NJ=NJ+1 

429  NPR I NT= 1 
MAS  =  1 

CALL  OUT  1 
L=NJ-1 

I F ( NC ( 14J.E0.0)  GO  TO  4295 

WR I TE { 26  *4  30 )  L.NJ.NJREG  * 

WRITE  TAPE  6  ALSO. 

4295  WR I TE ( 6*430 )  L  »NJ  *NJREG  * 

430  FORMAT { 33H0  REZi  IS  ABOUT  TO  COMBINE  Z0NES*I4.4H  AND»I4, 

1  10H  IN  REGION, 1 3  I 

COMBINE  ZONE(NJ)  WITH  ZONE(NJ-l>  NOTE  THAT  L+1=NJ 
WEIGHT  D »  E*  AND  P  ACCORDING  TO  MASS 

DMT=GRAMS1L)+GRAMS(L+1) 

X  <  NJ ) =X ( N J+ 1 J 

J=NJ 

NGEOM=2 

CALL  GEOM 
D <  L ) =DMT /VOL 

E(L)=(E(L)*GRAMS(L)+E(L+1)*GRAMS(L+1) J/DMT 
P(L)=(P(L) *GRAMS(L)+P(L+1 )*GRAMS(L+1 ) ) /DMT 
Q ( L )  =  ( 0 ( L ) +0 ( !  +1 i ) *  «  5 
GRAMS ( L ) *OMT 
PO ( L ) =P ( L ) +0 ( L  ) 

RELABEL  ZONES  BEYOND  ZONE  NJ. 

JLAST2=JLAST-2 
DO  440  L  =  N J » JLAST 2 
X( L)=X(L+1) 

U«L)=U(L+1 ) 

D  <  L ) aD ( L+  1 ) 

E ( L >=E(L+1) 

Q ( L ! =Q ( L+ 1  ) 

GRAMS ( L ) =GRAMS ( L  +  1  ) 

PO.DJDT ,DTZJ,TlMP,GAMMAJ  WILL  BE  RECALCULATED  ANYWAY  IN  MOTION 
AND  ARE  AC'  JSTEO  HERE  ONLY  TO  MAKE  THE  NEXT  PRINTOUT  CONSISTENT. 
PQ(U=PQ(L  +  1  ) 

DUDT(L)=DUDT  (L  +  l  ) 

DTZJ(L)=DTZJ(L+1) 
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TEMP ( L )* TEMP (L+l )  RE21 

GAMMAJ { L ) *GAMMAJ ( L+l } 

440  CONTINUE 

X< JLA5T1)*X{JLAST) 

U( JLAST1 )*U<  JLAST ) 

DUDT { JLAST1 ) *DUDT  C  JLAST } 


CALC  INTERJ  AND  SET  SOME  CONSTANTS 
NZONESC  NJREG ) *NZON£S(NJREG ) -1 
DO  445  I *NJREG ♦ IMAX 
INTERJf 1+1 )= INTERJ ( 1+1) -1 
445  CONTINUE 

JLAST=JLAST1 
JLAST1* JLAST2 
JCALC® JCALC-1 
JCALC1* JCALC-1 
NPRINT=1 
MAS*  1 

CALL  OUT  1 

I F ( KREZ3B-NJ )  499.490*490 
490  KREZ3B=KREZ38-1 
499  RETURN 
END 


B-38 


nrmonnnn 


NOL  TR  62-168 


$ I8FTC  REZ2  LIST *REF 
SUBROUTINE  REZ2 

THIS  SUBROUTINE  DOUBLES  THE  RANGE  OF  X.  WHILE  KEEPING  THE  NUMBER 
OF  ZONES  IN  THE  PROBLEM  THE  SAME.  PAIRS  OF  ADJACENT  ZONES 
ARE  AVERAGED  AND  REPLACED  BY  ONE  LARGER  ZONE# 

IF  REZ2  IS  USED*  THE  MAXIMUM  NUMBER  OF  ZONES  PERMITTED  IN  THE 

PROBLEM  IS  299  RATHER  THAN  300« 

FIRST  KREZ2-1  REGIONS  ARE  NOT  SUBJECT  TO  REZONING* 

*****  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 

C 

1FC IDENT.NE.l)  GO  TO  510 

502  FORMAT ( 40H0SUBROUT INE  REZ2  VERSION  3»  12-10-64  ) 

RETURN 

C 

510  JLAST4= JLAST— 4 

I F ( A8S(U( JLAST4J ) *LT*UCUT)  RETURN 

515  NPRINT=1 
MAS*1 

CALL  OUT1 

516  FORMAT?35HOREZONING  BY  REZ2  IS  ABOUT  TO  OCCUR) 

C 

L2  =  0 
1START=1 
NI X=KRFZ2-1 

c  DO  NOT  DO  ANYTHING  TO  THE  FIRST  NIX  REGIONS* 

IF(NIX)  521*521*518 
518  DO  520  I  =  1 » N I X 
L2=L2+NZONES< I ) 

520  CONTINUE 
LA*L2+1 
LB=L2+2 
ISTART=NI X+l 

C 

521  DO  550  I  =  1  ST  ART  *  I  MAX 
NZHALF=NZONES( I )/2 

IF (NZONESI I ) -2*NZHALF )  522*522*531 

5  on  anh  minT  WILL  BE  RECALCULATED  ANYWAY  IN  MOTION 

C  AND  ARE  DONE  HERE  ONLY  TO  MAKE  THE  NEXT  PRINTOUT  CONSISTENT. 

C 

522  IF(I-l)  599*523*524 

523  Ll=l 
L2=NZHALF 
LA*  1 
LB*2 

GO  TO  525 
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524  L1*L2+1 
L2=L2+NZHALF 
LA*LA+2 
L8*LA+1 

525  NZONESl I )*NZHALF 

C  COMBINE  ZONES  LA  AND  LB  TO  FORM  ZONE  L. 

X ( LI ) *X ( LA ) 

UCLl )*U(LA) 

DUDT (LI >*DUDT (LA) 

DO  530  L=L1 *L2 
DMT*GRAMS ( LA ) +6RAMS ( LB ) 

X ( L+l ) *X  <  LA+2 ) 

U(L+l)«U(LA+2) 

DUDT (LI) =DUDT { LA+2 } 

J*L+i 

NGEOM=2 

CALL  GEOM 
D( L ) «DMT /VOL 

E ( L  J  =  ( E ( LA  J  *GRAMS ( LA )  +E{ LB)*GRAMS(LB) ) /DMT 
P ( L )  =  ( P  <  LA ) *GRAMS ( LA }  +P ( LB ) *GRAMS ( Lb  > ) /DMT 
TEMP { L ) *— 0. 

DTZ J( L  >  =-0» 

Q( L ) * (0 ( LA ) +Q( LB  > )*.5 
GRAMS! L ) =DMT 
PQ ( L ) *P { L )+0( L  J 

GAMMA J { L ) « ( GAMMA J ( L A ) +GAMMAJ ( LB ) ) * . 5 
I F { L— L2 )  526,530.530 

526  LA*LA+2 
LBaLA+1 

530  CONTINUE 
GO  TO  550 

C  COMBINE  THE  FIRST  THREE  ZONES  OF  A  REGION  HAVING  AN  ODD  NUMBER 
C  OF  ZONES. 

C 

531  I F ( I  —  1 J  599  *535.537 

535  I F ( NZONES ( 1 >  — 1 >  599.548.536 

536  Ll=l 
L2  =  1 
LA  =  1 
LB=2 
LC  =  3 

GO  TO  538 

537  IF(NZONES( I )-l )  599,5380.5370 
5380  LI *L2+1 

L2=L2+1 
LA*LA+1 
LB*LB+1 
GO  TO  550 
5370  LI =L2+1 


REZ2 


B-40 


o  n  o  n 


NQL  tr  62-168 


L2*L1 
LA*LB+1 
LB*LA+1 
LC*LA+2 

538  DMT*GRAMS ( LA ) +GRAMS ( LB ) +GRAMS ( LC ) 

X  t  Ll  +  1 ) *X ( LA+3 ) 

U( Ll+1 ) *U( LA+3 ) 

J=L1+1 
NGE0M*2 

CALL  GEOM 
D( LI )*DMT/VOL 
U( LI ) *U( LA) 

DUDT { LI ) *DUDT ( LA ) 

Q( LI )  =  (QILA )+Q ( LB ) +Qt  LC ) ) /3» 

E(L1 )*{E(LA)*6RAMS(LA)  +£ ( LB ) *GR AMS  C  LB )  +E(LC)*GRAMSCLC) l/DMT 
P ( LI ) = ( P ( LA)*GRAMS( LA )  +P ( LB J+GRAMS < LB )  +P(LC)*GRAMStLCi >/DMT 
TEMP(Ll>*-0. 

DTZJ(L1)=~0. 

GRAMSIL1 )*DMT 
PQ(L1)=P(L1)+Q(L1) 

GAMMAJ <  L 1 )  =  < GAMMA J { LA ) +GAMMAJ <  LB ) +GAMMAJ ( LC ) ) /3 • 

NZONES! I ) =NZHALF 
I F ( NZONES ( I ) — 1 )  541,541*540 

540  L1=L1+1 
L2*L1+NZHALF— 2 
LA*LC+1 
LB*LC+2 
GO  TO  525 

541  L2=L1 
LA«LC-1 
GO  TO  550 

548  Ll*l 
L2*l 
LA=0 
LB*  1 

550  CONTINUE 

CREATE  A  NEW  REGION  UMAX+l).  NEW  REGION  BEGINS  WITH  ZONE  LNEW 
ITS  INNER  BOUNDARY  IS  INTERFACE  LNEW. 

LNEW=L2+1 

NZONES ( IMAX+1 ) *JLAST-LNEW 
XCONST=X ( LNEW-1 )-X(LNEW-2) 

C  USE  ARDC  AMBIENT  CONDITIONS  FOR  NEW  ZONES  IF  ARDC  IS  USED. 

JLAST1=JLAST-1 
X ( LNEW ) *X ( LNEW-1 )+XCONST 
1F( IARDC.GE.31 )  GO  TO  561 
555  DO  560  L=LNEW* JLAST 1 
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X(L+1 ) =X(L)+XCONST  REZ2 

ALTCM«COSPHI*(X(L)+X(L+m  +HOBKM*l*E+5 
CALL  AROC(ALTCM.PRATIO»TRAT10»ORATIO) 

D(L)*DRATIO*. 001225 
E  ( L ) »PRAT 10*1 #01 325E6/D  <  L; /,4 

560  CONTINUE 
GO  TO  563 

ZONE  JLAST1  STILL  CONTAINS  AMBIENT  DATA  FROM  BEFORE  REZONING, 

561  DO  562  L*LNEW» JLAST1 
D(L)=D( JLAST1) 

E(L)*E( JLAST1) 

X ( L ) -X ( L-l >+XCONST 

562  CONTINUE 

563  DO  564 ' L=LNEW» JLAST 1 
PQ(L)=PQf JLAST!) 

TEMP ( L ) *-Q# 

GAMMAJI L ) =-0« 

DTZJ ( L ) *-0« 

Q(L)«0. 

U(L)=U( JLAST) 

564  CONTINUE 

X ( JLAST ) =  X ( JLAST 1 ) +XC0NST 

571  DO  572  L*LNEW* JLAST1 
NGE0M=2 

J=L+1 

CALL  GEOM  - _ _ 

GRAMS(L)xD(L)*VOL 

572  CONTINUE 

NEQST 1  I MAX+1 ) *NEQST  C IMAX  J 
C I  NQ  ( IMAX+1 )  *=C  I NQ  (  IMAX  ) 

A I  NO ( j  M4X+ 1 ) * A I  NO { IMAX) 

GAMMA I ( IMAX+1 ) “GAMMA  I (IMAX) 

IMAX*IMAX+1 

OUTBDY ( IMAX ) *X ( JLAST ! 

JUTEST= i JUTEST-LZONES ) /2+L ZONES 
DO  573  I *1 ♦ IMAX 

573  INTERJ( 1+1 )= INTER J( I )+NZONES( I ) 

TLISTt 1 )=2.5*TLIST<  1) 

TLIST(3)=2.5*TLIST(3) 

TLIST(5)=2.5*TLIST( 5) 

TLIST(6)=2«5*TLIST(6) 

NPRINT=1 

MAS  =  1 

CALL  OUT  1 
599  RETURN 
END 
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$  I BFTC  REZ3  LIST .REF 

SUBROUTINE  REZ3 

THIS  SUBROUTINE  MAINTAINS  FINE  ZONES  IN  THE  SHOCK  FRONT* 
*#*##CAUTION*****  REZ3  CANNOT  BE  USED  IN  SAME  PROBLEM  WITH  REZ2. 

**»  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 

EQUIVALENCE  I  NELL .KREZ3 ) * ( NANCY  » KREZ3A ) • (NOLA  .KREZ33) 

NELL  =  NUMBER  OF  FINE  ZONES  PER  LARGE  ZONE. 

NANCY  =TOT AL  NUMBER  OF  FINE  ZONES. 

NOLA  =OUTER  INTERFACE  OF  OUTERMOST  FINE  ZONE. 

NOREEN=  NUMBER  OF  ZONES  BACK  FROM  NOLA  FOR  MAKING  REZONE  TEST. 

IF( IDENT.NE.l)  60  TO  603 

601  WR I TE ( 6*602 ) 

602  FORMAT ( 55H0SU8ROUT INE  REZ3  VERSION  3.  3-  1-65  TAPES  6*  26 

MZOAN=0 
NOREEN=3 
RETURN 

603  NORA=NOLA-NOREEN 

TRIGGER  REZONING  IF  ZONE  OF  MAX  Q  IS  AT  OR  BEYOND  ZONE  NORA. 
FIND  WHERE  QMAX  IS. 

JMAXQ= 1 
QMAX=0 • 

DO  605  J=1 ♦ JCALC 
I F ( Q C J) .LE.QMAX)  GO  TO  605 
60*  QMAX=Q( J } 

JMAXQ= J 
605  CONTINUE 

I F ( JMAXQ.GE.NORA )  GO  TO  607 

IF  NC ( 8 ) =0 ♦  TRIGGER  REZONING  IF  PQ(NORA)  EXCEEDS  PQ<NORA+l) 

BY  20  PERCENT. 

IF(NC(8).NE.O)  RETURN 

IF(PQINORA) „LT.1.2*PQ(NORA+l) )  RETURN 

607  JST ART  =  NOLA-NANCY 
JEND=JSTART+NELL~1 

C  PRINT  BEFORE  AND  AFTER  EACH  REZONING  IF  NC(16)*1. 

I F ( NC ( 16 ) .NE.l )  GO  TO  611 

608  NPR I NT= 1 
MAS=  1 

CALL  OUT  1 
C 

C  FORM  ONE  LARGE  ZONE  FROM  FIRST  NELL  FINE  ZONES 

611  ZNELL=NELL 
PTOT=0. 

GTOTxO. 

T0TIE=0» 
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70TKE*0. 

615  DO  620  J*JSTARY.JEND 
PTOT*PTOT+P( J) 

6T0T  *GTOT+GRAMS  f  J ) 

USG* .5*  (U(  J) **2+U ( J+). ) **2  ) 

TOT  1 E*TOT I E+E ( J ) *GRAMS ( J ) 
rOTKE*TOTKE+. 5«USQ*GRAMS ( J ) 

620  CONTINUE 

TOTE*TOTiOTOTlCE 

GRAMS(JSTART)eGTOT 

X(JSTART+1)*X(JEND+1) 

J«JSTART+1 

NGE0M*2 

CALL  6E0M 

D( JSTART )«GR AMS {JSTART ) /VOL 
P{ JSTART )*PTOT /ZNELL 
0( JSTART)*0. 

'  PQ ( JSTART )  =P (JSTART ) 

C  ADJUST  INTERNAL  ENERGY  SO  AS  TO  CONSERVE  TOTAL  ENERGY. 

JPNELL* JSTART+NELL 

EKJST  *  ( U( JSTART ) **2+U{ JPNELL ) **2 )*GRAMS ( JSTART )  /  4. 

El JST*TOTE~EKJST 
E  <  JSTART  )*EUST/GRAMS<  JSTART) 

C 

C  RELABEL  EXISTING  FINE  ZONES 

JENDP= JEND+1 
N0LA1*N0LA— 1 

634  DO  635  J= JENDP *N0LA1 
JNEW* J-NELL+1 

X ( JNEW ) *X { J ) 

U( JNEW ) *U<  J ) 

D( JNEW)=D( J) 

P( JNEW)*P( J) 

Q( JNEW) =Q( J ) 

E( JNEW)=E( J) 

GRAMS! JNEW) *GRAMS ( J ) 

PCM  JNEW) *PQ( J) 

635  CONTINUE 
JNEW-NOLA-NELL+l 
X( JNEW)=X(NOLA) 

U< JNEW) *U< NOLA) 

C 

C  DIVIDE  NEXT  LARGE  ZONE  INTO  NELL  FINE  ZONES. 

C  NOLA  STILL  REFERS  TO  THE  LARGE  ZONE  ABOUT  TO  BE  DIVIDED. 

640  JDONE*NOLA+l 
UAMB*U( NOLA+1 ) 

DAMB*D< NOl  A ) 

PAMB*P(NOLA) 

EAMB*E( NOLA) 

XS I ZE*(X< NOLA+1 )-X< NOLA) ) /ZNELL 
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C  ADJUST  INTERPOLATION  IF  NOLA+1  IS  LAST  INTERFACE  IN  PROBLEM. 

IF(NOLA.GE.JLAST-l)  60  TO  642 

641  NYMPH*N0LA+1 
GO  TO  643 

642  NYMPH=NOLA~  1 

643  XI A= ( X ( NOLA ) +X ( NOLA+1 ) ) /2 . 

XI ( X ( NYMPH ) +X ( NYMPH+1 )} /2 • 

X1D*X1B-X1A 

C 

JBEGIN*N0LA-N£LL+2 
DO  660  J=JBEGIN.JDONE 
J=J 

JMH* J— 1 

X  <  J ) *X ( J-l I+XSIZE 
U< J)*UAM8 

I F ( I ARDC.LT .31 )  GO  TO  6431 

6430  P«JMH)*PAMB 
D< JMH)=DAM8 
E( JMH)*EAM3 
GO  TO  6450 

C  INTERPOLATE  WITHIN  THE  LARGE  20NE  IF  ARDC  IS  USED. 

6431  DINT=D(NYMPH) 

PINT*P ( NYMPH) 

EINT*E( NYMPH) 

XlC=X(J)-X5IZE/2. 

XI R*  ( X1C— X1A ) /X1D 

645  D( JMH )=DAMB+ ( DINT-DAMS )  *X1R 
P< JMH)=PAMB+(PINT-PAMB)*X1R 
E  <  JMH )  =  EAMB+ ( E I NT-EAMB )  *X  1R 

6450  Q( JMH ) =0. 

PQ ( JMH) *P( JMH ) 

PS  I ( JMH ) *PQ ( JMH  >  *14. 50382E-6 
TEMPI JMH) *0. 

DTZJ  < JMH) =DTZJ ( NOLA ) 

GAMMA J ( JMH ) *GAMMAJ ( NOLA ) 

C 

646  NGE0M=2 

CALL  GEOM 

GRAMS ( JMH ) *D ( JMH ) * VOL 
660  CONTINUE 

ADJUST  NZONES  AND  INTERJ. 

1000  NELL1*NELL-1 

C  NLAST  IS  INNER  BOUNDARY  OF  FINE  ZONES. 

NLAST«NOLA-NANCY 
DO  1020  1 1 *2  *  I  MAX 
IFINOLA.LT .INTERJ (II))  GO  TO  1020 
IF (NLAST .GE* INTERJ (II))  GO  TO  1020 
1002  INTERJ! 1 1 )»INTERJ( 1 1 ) -NELLI 
NZONES ( 1 1 ) =NZ0NES ( 1 1 J+NELL1 
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NZONES! 1 1-1 }«NZONES< 1 1-1 I-NELL1  REZ3 

1020  CONTINUE 
C 


661  NOLA-NOLA+I 
MZ0AN«MZ0AN*1 
IF ( NC! 14 ) *EQ«0  >  GO  TO  665 

WRITE! 26*666 )  MZOAN.NOLA .NCYCLE  * 

GO  TO  6660 

665  WRITEC6.666)  MZOAN.NOLA  »NCY<  * 

666  FORMAT! 25N0REZ0NING  BY  REZ3»  M2wrtN=I3.6H  NOLA»I3*8H  NCYCLE*I4) 

6660  IF INOLA-JLAST )  669.667.667 

667  KREZ3*0 

IF!NC!14).EQ.0>  GO  TO  6670 

WRITE! 26.668  I  * 

C  WRITE  TAPE  6  ALSO. 

6670  WRI TE ( 6.668 )  * 

668  FORMAT! 36H1FINE  ZONES  HAVE  REACHED  OUTER  WALL.  /53H0  KREZ3  IS  NOW 
1  ZERO  TO  PREVENT  FURTHER  USE  OF  REZ3.  ) 


669  IF(NC! 16I.NE.1 )  RETURN 

670  NPRINT*1 
MAS*  1 

CALL  OUT  1 
699  RETURN 
END 
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S1BFTC  REZ4  LIST .REF 

SUBROUTINE  REZ4  REZ4 

REZONE  FIRST  REGION 

THIS  SUBROUTINE  REDUCES  THE  NUMBER  OF  ZONES  IN  REGION  1  TO  XREZ4 
ZONES  AT  TIME  TREZO.  REGION  1  OFTEN  KEEPS  THE  TIMESTEP  SMALL  AFTER 
ALL  PHENOMENA  OF  INTEREST  HAVE  MOVED  OUT  INTO  OTHER  REGIONS*  THIS 
ROUTINE  PERMITS  THE  PROBLEM  TO  RUN  FASTER  IN  SUCH  CASES. 

***  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 

IFC IDENT.NE.l)  GO  TO  403 

401  WRITE! 6*402 )  * 

402  FORMAT ( 40H0SUBROUT INE  REZ4  VERSION  3,  3-  1-63  ) 

ND0NE=0 
RETURN 

DO  NOT  PROCEED  UNTIL  T  HAS  REACHED  TREZO. 

403  IFIT.LT. TREZO)  RETURN 

C  DO  NOT  PROCEED  UNLESS  PROBLEM  HAS  SPHERICAL  SYMMETRY. 

I F ( K.NE.3 )  RETURN 

C  DO  NOT  PROCEED  IF  KREZ4  IS  NOT  LESS  THAN  NZONESCl). 

IF(NZ0NES(1).LE.KREZ4)  RETURN 
405  EKIN1=0» 

PTOT=0. 

GTOT=0. 

TOT  I E=0. 

TOTKE=0. 

J5TOP=NZONES ! 1 ) 

ZN=NZ0NES(1) 

NPR INT=1 
MAS=  1 

CALL  OUT  1 

WRITE(6»415)  JST0P.KREZ4.NCYCLE  * 

415  FORMAT! 37H0REG ION  1  IS  ABOUT  TO  BE  REDUCED  FROM,I4»3H  TO.I4.15H  ZO 
1NES  AT  CYCLE. 15) 

DO  420  J= 1  * JSTOP 
PTOT  =PTOT  +  P ( J) 

GTOT=GTOT+GRAMS! J) 

USQ=.5*(U( J)**2+U! J+l)**2) 

TOT  I E* TOT  I E+E ( J ) *GRAMS ( J  > 

TOTKE= TOTKE+ . 5*USQ*GRAMS ( J ) 

420  CONTINUE 
C 

ETOTl=TOTIE+TOTKE 
1 1 1  =  INTERJ ! 2 ) 

XXX  =  X ( I  I  I ) 

ZZZ=KREZ4 
SI ZE*XXX/ZZZ 
V0L1»4. 1887902*XXX**3 
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DAV=GTOT /V0L1  REZ4 

PAV*PTOT/ZN 
00  430  J=1  *KREZ4 
X< J+l )=X( J)+SIZE 

VOLUMES,  1887902*  < X < J+l ) **3-X < J ) **3  ) 

GRAMS ( J ) *0 A V* VOLUME 
P< J)=PAV 
Q< J)=0. 

UC J)=U< JSTOP)*X< J)/XXX 
D  (  J )  =DAV 
DTZJ( J)=-0. 

TEMP ( J ) =-0« 

PO(J)=P( J) 

EKIN1=EKIN1+(U( J)**2+U( J+l )**2>*GRAMS< J)/4. 

430  CONTINUE 

ADJUST  E  TO  CONSERVE  TOTAL  ENERGY  IN  REGION  1  • 

EINT1=ET0T1-EKIN1 
EP£RG=E INTl/GTOT 
DO  440  J=1*KREZ4 
440  E ( J } =EPERG 

RELABEL  REST  OF  ZONES 
JSTART=KREZ4+1 
JC=JST0P-XREZ4 
JEND=JLAST-JC 
NZONESI 1)=KREZ4 
DO  450  J=JSTART » JEND 
JW= J+JC 
X< J)*X< JW) 

U( J)=U( JW) 

D( J! =  D( JW) 

E( J)=E( JW) 

Q< J)=Q( JW) 

PQ( J)=PQ{ JW) 

GRAMS'’  J  )  =GRAMS  ( JW  ) 

GAMM>'  (,  J)=GAMMAJ(  JW) 

DTZJ( J)=DTZJ( JW) 

DUDT ( J ) =  DUDT ( JW ) 

DWAS ( J ) =DWAS ( JW ) 

EWAS( J)=EWASJ JW) 

PWAS ( J ) =PWAS  <  J W ) 

TEMP ( J ) *TEMP  < JW ) 

450  CONTINUE 
C 

IMAX 1= I MAX+1 
DO  460  1=2* IMAX1 
460  INTERJf I )=!NTERJ( I )-JC 
JLAST*=JLAST-JC 
JLAST1=JLAST-1 
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JCALC*JCALC-JC 
KREZ3B=KREZ3B~JC 
NPRINT=1 
MAS*  1 

CALL  OUT 1 

C  PREVENT  FURTHER  USE  OF  THIS  ROUTINE  BY  SETTING  KREZ4«0* 

KR£Z4=0 
499  RETURN 
END 
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SIBFTC  ARDC  LIST  »REF 

SUBROUTINE  ARDC ( ALTCM.PRAT IO» TRAT 10, DRAT  10 ) 

REFERENCE 

MINZNER*  R.A.  AND  RIPLEY*  W.S.*  THE  ARDC  MODEL  ATMOSPHERE* 
1956*  AFCRC  TN-56-204*  AD-110233*  DEC  56*  UNCLASSIFIED. 
ALTCM*ALT I TUDE  IN  CENTIMETERS 

PRAT IO=RAT 10  OF  AMBIENT  PRtSSURE  TO  SEA  LEVEL  PRESSURE  (1  ATM) 

TRAT IO=RAT 10  OF  AMBIENT  TEMPERATURE  TO  SEA  LEVEL  TEMP  (288.16 
DRAT 10=RATI0  OF  AMBIENT  DENSITY  TO  SEA  LEVEL  DENSITY! .001225  6/CC) 

100  ALTZ»ALTCM/100. 

ALTH=6356766.0*ALTZ/( 6 356766. Q+ALTZ) 

IF(ALTH. 6T. 11000. )  GO  TO  102 

101  TEMP=288.16-0.0065*ALTH 

PAMB*14. 6961 78/ (288.1 60/ ( 288. 160-0. 0065*ALTH) ) **5.25612218 
GO  TO  118 

102  I F ( ALTH.GT .25000. )  GO  TO  104 

103  TEMP=216.66 

PAMB=3. 28254528/ < 10. ** < 0.068483253E-3* ( ALTH-1 1000. OH ) 

GO  TO  118 

104  IF ( ALTH.GT .47000. )  GO  TO  106 

105  TEMP=216«66+0.003* ( ALTH-25000.0 ) 

PAMB=0. 36094654/1  ( 141 .660+3. 0E-3*AL.TH  > /2 16.66 )  **1 1 . 38826473 
GO  TO  118 

106  IF( ALTH.GT. 53000. )  GO  TO  108 

107  TEMP=282.66 

PAMB=0. 0174686/(10. **(0. 0524926823E-3*(ALTH-47000. 0) ) ) 

GO  TO  118 

108  IF( ALTH.GT. 79000*)  GO  TO  110 

109  TEMP=282. 66-0. 0045* (ALTH-53000.0) 

PAMB=8.40408E-3/( ( 282. 66/TEMP ) **7 . 592176 ) 

GO  TO  118 

110  IKALIH.OI  .90000. )  GO  TO  112 

111  TEMP=165.66 

PAMB=1.46198E-4*EXP  (-0.0341647942* ( ALTH-79000.0) / 165. 66 ) 

GO  TO  118 

112  IF(ALTH.GT. 105000.)  GO  TO  114 

113  TEMP=165.66+0.0040»( AL1H-90000.0)  '  ' 

PAMB= 1.5519E-5*( 165 .66 /TEMP )**8.541198 
GO  TO  118 

114  IK ALTH.GT. 160000. )  GO  TO  116 

115  TEMP =22 5 .66+0. 02 * ( ALTH-105000.0 ) 

PAMB=1.04442t-6^( 225. 66/TEMP) **1.708239 
GO  TO  118 

116  IF(ALTH. GT. 170000. )  GO  TO  119 

117  TEMP= 1325. 66+0.01* < ALTH-160000.0) 

PAMB=5o 1401 5E-8*( 1325 066/ TEMP >**3.4164794 
GO  (0  118 

liy  IK  ALTH.GT  .200000.  )  GO  TO  121 
120  TEMP* 1425 • 66+0.005*(ALTH- 170000.0) 
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PAMb=4*0664t-8#( 1426* 66/ TtMP)##6* 832958 
GO  TO  118 

121  TEMP= 1575 *66+0* 0035# (ALTH-2OQOO0*O ) 

PAMB=2*0595E-8# ( 15  76* 66 /TEMP ) ##9* 76 1369 
118  TRATI0=TEMP/288.16 

PRATI0*PAMB/14. 696178 
DtNSAM=3« 2 365427  fc.-4*PAMB/ TEMP 

C  OtNSAM  IS  IN  UNITS  OF  SLUGS  PER  SQUARE  INCH  PER  FOOT* 
DRAT IO=DENSAM/ *0000 165 
RETURN 
END 
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SIBFTC  GtOM  LIST *REF 

SUBROUTINE  GtOM 
C 

*****  INShKI  SIANUAKU  COMMON  CAKOS  btFORc  COMPILING  ***** 
C 

IM  IDtNT.Nt.l)  GO  TO  9 

5  WR  I.TE  (6*6) 

6  FORMAT (40H0SUBRGUTINE  GEOM  VERSION  3,  3-  1-65  J 

RETURN 


CALC  AREA  OF  INTERFACE  J  OR  VOLUME  OF  ZONE  J-l 
9  NGEOM*NGEQM 

GO  TO  ( 10*20) *NG£OM 

10  K«(C 

GO  TO  < 11 »12»13)  »K 

11  AREA«1. 

GO  TO  99 

12  AREA*6.2831853*X(J) 

GO  TO  99 

13  AREA* 12 • 566372 *X <  J ) **2 
GO  TO  99 

20  K*K 

GO  TO  ( 21 *22  *23 ) *K 

21  VOL*X( J)-X( J-l ) 

GO  TO  99 

22  V0L*3. 1415927* (X( J)**2~X< J-l)**2) 

GO  TO  99 

23  VQL*4. 1887902* <X(J)**3-X< J-l 5**3) 

99  RETURN 

END 
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SI8FTC  EQST  LIST, REF 
SUBROU I  I  ME  EOS  I 
C 

*****  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 

AN  EQUATION  OF  STATE  MUST  CALCULATE  P(J)»  E<J),  AND  SOUND(J) 
AND,  IF  DESIRED,  TEMP^J)  AND  GAMMA J ( J )  • 

NQST*NQST 

GO  TO  (1,2, 3, 4, 5,6, 7,8, 9,10, 11, 12) ,NQST 

1  CALL  EQS1 
GO  TO  100 

2  CALL  EQS2 
GO  TO  100 

3  CALL  EQS3 

GO  TO  100 

4  CALL  EQS4 

GO  TO  100 

5  CALL  EQS5 

GO  TO  100 

6  CALL  EQS6 

GO  TO  100 

7  CALL  EQS7 

GO  TO  100 

6  CALL  EQS8 

GO  TO  100 
9  CALL  EQS9 

GO  TO  100 

10  CALL  E010 

GO  TO  100 

11  CALL  EOl 1 

GO  TO  100 

12  CALL  EQ12 

100  RETURN 

END 
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S18FTC  EQS1  LIST. REF 

SUBROUTINE  EQS1  EOS1 

EQUATION  OF  STATE  FOR  GAMMA-LAW  GAS.  E»PV/!G-1>. 

#**  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 

AN  EQUATION  OF  STATE  MUST  CALCULATE  P.  E.  AND  SOUND. 

IF  DESIRED.  IT  MAY  ALSO  CALCULATE  TEMP  AND  GAMMAJ. 

IF C IDENT.NE.l)  GO  TO  800 
IF!NOMOR£.EQ.71)  RETURN 

101  WRITE ! 6. 102 )  * 

102  FORMAT ! 4CHOSU8ROUT I NE  EQS1  VERSION  3,  3-  1-65  ) 

RGAS«2.87043E.v6 
N0M0Rt*71 

103  RETURN 
C 

800  GAMMA*GAMMAI ( I ) 

C  BEGIN  ENERGY  ITERATION. 

90  GA*GAMMA-1 .0 

GB* ( GAMMA+1 • } /GA 
Vl*l./D! JMH) 

V2«1./DWAS!JMH) 

VB*V1— V2 
VC»GB*V1-V2 

P < JMH >  = « 2 . *EWAS <  JMH )  - ! PWAS < JMH ) +2 • *Q ( JMH )  ) * VB ) / VC 
PV«P! JMH5/D! JMH) 

E ! JMH ) *PV/GA 
TEMP ( JMH ) =PV/RGAS 
GPV*GAMMA*PV 
IF(GPV)  211.211.220 

211  WRI TE (6.212 )  JMH . E ! JMH ) »D ( JMH ) .P ( JMH ) .GAMMA 

212  FORMAT ( 18H1*****ST0P  IN  EQS1/6H  Z0NE=I4/8H  E( JMH) «1PE14.5/8H  D( JMH 
1 ) *1PE14.5/8H  P( JMH)=1PE14.5/8H  GAMMA  -1PE14.5) 

KUT0FF=1 
MAS  =  1 

CALL  OUT  1 

220  SOUND! JMH) *SQRT (GPV) 

GAMMAJ! JMH)*GAMMA 

RETURN 

END 
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SIBFTC  EQS2  LIST »REF 

SUBROUTINE  EQS2  E0S2 

APPROXIMATE  EQUATION  OF  STATE  FOR  REAL  AIR  (VARIABLE  GAMMA) • 

***  INSERT  STANDARD  COMMON  CARDS  BEFORE  COMPILING  ***** 

DATA  RGAS/2.87043E+6/ 

IF( IDENT.NE.l)  GO  TO  800 
I F ( NOMORE. EQ. 73 )  RETURN 

101  WRITE! 6.102 )  * 

102  FORMAT (45H0SUBR0UTINE  EQS2  VERSION  3,  3-  1-65  AIR  ) 

NOMORE=73 

103  RETURN 

CHANGE  ENERGY  FROM  ERG/GRAM  TO  MBAR  CC/GRAM. 

800  ENERGY=E(JMHS*1.E-12 
DENRAT=D( JMH)/. 001225 
I F ( ENERGY .GT.O. )  GO  TO  4 

202  WRI TE( 6*212 )  JMH*E( JMH) »D( JMH) »P{ JMH) *GAMMA 

212  FORMAT ( 18H1*****STQP  IN  EQS2/6H  Z0NE=I4/8H  E( JMH) *1PE14.5/8H  D{ JMH 
1 ) -1PE14«5/8H  P(JMH)=1PE14.5/8H  GAMMA  *1PE14.5) 

KUT0FF*1 
MAS-1 

CALL  OUT  1 

STOP 

4  IF(ENERGY.GT. 0.003)  GO  TO  204 

2  GAMMA* 1.4 

3  GO  TO  90 

204  ELOG=ALOG10 ( ENERGY ) 

IF ( DENRAT )  202.202.205 

205  DLQG*ALOG 10 (DENRAT ) 

6  IF(ENERGY.GT. 0.015)  GO  TO  8 

7  GAMMA*-. 071533*ELOG+1. 21953 
GO  TO  90 

8  IF ( ENERGY. GT .0*05 )  GO  TO  10 

9  AG*-. 172124 
BG=1. 03606 
CG*-. 338509 
DG=. 732590 
GO  TO  80 

10  IF ( ENERGY.GT .0.3 )  GO  TO  12 

11  AG*-. 109233 
BG*1. 11789 
CG*-. 106663 
DG*1. 03423 
GO  TO  80 

12  IF( ENERGY.GT .20* )  GO  TO  14 

13  AG*. 0411205 
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8G* 1.1965  £052 

CG-.0137068 
DG* 1.09 72 
GO  TO  80 

14  IF  $  ENERGY. GT» 70. )  GO  TO  16 

15  AG*. 27569 
8G*. 89132 
CG*.52383 
DG*. 43348 
GO  TO  80 

16  GAMMA*. 086588*ELOG+ 1.24024 
IF(GAMMA.GT« 1.667}  GAMMA*1,667 
GO  TO  90 

80  GAMI*AG«ELOG+BG 
GAM2*CG*£L0G+0G 

GAMHA*GAM1- ( l.-DLQG) /7.*( GAM1-GAM2 ) 

CALCULATION  OF  AN  EFFECTIVE  GAMMA  NOW  COMPLETED. 

BEGIN  ENERGY  ITERATION. 

90  GA*6AMMA-1 .0 

GB*(GAMMA-H.)/GA 
Vl*l./D( JMH) 

V2*1./DWAS( JMH) 

DV*V1-V2 

VC«GB*V1-V2 

P<  JMH)*(2.*EWAS(JMH)~(PWAS«  JMH)+2.*QUMH)  )«DV)/VC 

22  PV*P<JMH)/D< JMH) 

E( JMH)=PV/GA 
TEMP ( JMH) *PV/RGAS 
GPV*GAMMA*PV 
IF(GPV)  202.202*221 
221  SOUND! JMH) * SORT (GPV ) 

GAMMAJ ( JMH ) *GAMMA 
999  RETURN 
END 
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